Atelier R ANF FOSS : La régression géographiquement Pondérée ou GWR
Objectif de l’atelier : Être en mesure de réaliser et interpréter une régression géographiquement pondérée (GWR).
Prérequis : Connaissances de base en analyse de données statistiques, avec au moins une compréhension de ce qu’est une corrélation et une régression.
1 Pourquoi la GWR ?
Lorsque l’on souhaite dépasser la simple caractérisation d’attributs liés à des individus statistiques, on fait appel à des méthodes de modélisation explicitant des relations statistiques. La méthode la plus employée, principalement en sciences sociales, pour mesurer et analyser la nature des relations et des effets entre deux ou plusieurs variables, est le modèle de régression linéaire.
Le principe de la régression linéaire est de modéliser la variable que nous souhaitons étudier (aussi appeler variable dépendante, VD) comme une fonction linéaire des variables que nous aurons définies comme explicatives de la VD (aussi appelées variables indépendantes, VI). Lorsque l’on s’intéresse à un phénomène social avec une emprise sur un espace, la régression linéaire pose plusieurs problèmes :
Le premier est empirique. La régression linéaire nous permet d’obtenir des coefficients (appelés betas β) et des résidus (notés epsilon ε). Ces β représentent l’effet de nos VI sur notre VD. Ces β sont considérés comme globaux, sans variation. Autrement dit, les modèles de régression linéaires considèrent que les VI interviennent de la même manière et avec la même importance sur l’ensemble de notre jeu de données. Si cette hypothèse peut être validée sur des populations statistiques définies aléatoirement et sans effet de structure a priori des VI ou de la VD, elle n’est que rarement vérifiée sur des données spatiales. En effet, les caractéristiques propres de chaque territoire (l’unicité de chaque lieu) impliquent que l’effet constaté en un lieu n’est pas forcément valable en un autre lieu de l’espace.
Concernant les prix des valeurs foncières que nous détaillerons par la suite dans cette fiche, on peut comprendre que la proximité au littoral, très prégnante en certains points de l’espace, ne joue absolument aucun rôle dans d’autres lieux. De même, certaines caractérisations du monde rural n’interviennent plus lorsqu’on se situe dans des milieux fortement urbanisés. Ainsi, les données spatialisées sont soumises à l’hétérogénéité spatiale : l’effet de nos VI va varier en fonction de l’espace. Un coefficient qui serait global et uniforme pour mesurer un effet paraît plus simple et donc tentant, mais non pertinent en géographie ; sur ce point nous pouvons nous référer à l’article de Brunsdon, Fotheringham et Charlton (Brunsdon et al. 1996). Ce concept d’hétérogénéité dans l’espace se traduit en statistique par celui de non stationnarité.
Le deuxième problème est statistique : chaque méthode statistique doit répondre à un certain nombre de conditions de validité. La régression linéaire ne fait pas exception. Trois conditions doivent être validées pour qu’une régression linéaire puisse être effectuée sans que l’interprétation des résultats conduisent à des raisonnements fallacieux :
- Les individus statistiques doivent être indépendants
- Les résidus doivent suivre une distribution normale
- Il ne peut pas y avoir plus de VI que d’individus statistiques
Si les deux dernières conditions ne trouvent pas de matérialisation spécifique sur des données spatiales, la première quant à elle concrétise un problème récurrent sur les données en géographie. Par leur nature même, les données spatiales ne peuvent pas remplir cette condition fondamentale pour une régression classique. La première loi de la géographie de Tobler : “everything is related to everything else, but near things are more related than distant things” en est une traduction tout à fait parlante.
Le troisième problème est lié aux effets de contexte : on ne peut pas étudier des données spatiales sans considérer que les individus statistiques (les objets spatiaux) appartiennent eux-mêmes à des agrégats qui ont une influence sur la variable à expliquer. Ainsi, certains attributs de l’agrégat vont avoir une influence sur l’entité spatiale de cet agrégat.
Le quatrième problème est lié à la problématique du MAUP (Modifiable Area Unit Problem). Il concerne un problème d’échelle d’application de la régression et peut conduire à des observations erronées. Une corrélation constatée à une échelle peut être uniquement liée à l’agrégation réalisée à cette échelle, mais s’avérer erronée à une échelle plus fine, invalidant de fait la relation entre les phénomènes étudiés (Mathian, Piron 2001). Certaines agrégations peuvent également varier et cacher des relations entre individus (Bailey, Gatrell 1995).
La GWR ne répond pas à l’ensemble de ces problèmes mais va nous permettre de résoudre les deux premiers en intégrant la dimension spatiale de nos données tout en tenant compte de l’hétérogénéité (ou non stationnarité) de leur effet.
1.1 Les packages
Voici les packages que nous utiliserons :
# Chargement, visualisation et manipulation de la données
library(here)
library(DT)
library(dplyr)
# Analyse et représentation statistique
library(car)
library(correlation)
library(corrplot)
library(ggplot2)
library(gtsummary)
library(GGally)
library(plotly)
# Manipulation et représentation de la données spatiales (cartographie)
library(sf)
library(mapsf)
library(rgeoda) #permet en plus de calculer les indices d'auto-corrélation spatiale
library(RColorBrewer)
# Calcul du voisinage et réalisation de la GWR
library(spdep)
library(GWmodel)1.2 Présentation et préparation des données
Nous avons cherché à traiter ici une variable présentant des caractéristiques spatiales fortes et qui rencontrent les deux problèmes exposés précédemment d’indépendance statistique et de non-stationnarité. Les prix de l’immobilier en France sont effectivement soumis à ces problèmes et peuvent être expliqués, au moins partiellement, par des variables quantitatives issues de données de l’INSEE. Nous avons traité l’information liée au prix de l’immobilier à l’échelle de Etablissements Publics de Coopération Intercommunale (EPCI) pour nous assurer un nombre d’observation suffisant dans chaque entité spatiale.
- Les données du prix de l’immobilier par EPCI (prix médian au m²) sont issues des ventes observées sur l’année 2018, extraites depuis la base de données des notaires de France par Frédéric Audard et Alice Ferrari. Ce fichier a été simplifié pour ne conserver que les variables d’intérêts parmi une cinquantaine
- Les données statistiques proviennent de l’INSEE (année 2019) : 9 variables ont été choisies pour leur potentialité à expliquer les variations des prix de l’immobilier, concernant la population, le logement et les revenus et niveaux de vie. Elles sont détaillées un peu plus bas.
- Comme indiqué en introduction les données spatiales proviennent de la base ADMIN-EXPRESS de l’IGN en accès libre. La couche est utilisée est celle des EPCI de la base ADMIN-EXPRESS-COG édition 2022 par territoire pour la France métropolitaine.
Les données statistiques et du prix de l’immobilier ont été
regroupées dans un même fichier CSV donnees_standr.csv.
Les données spatiales (géométrie des EPCI) sont au format shapefile dans le
fichier EPCI.shp.
1.2.1 Chargement des données sur le prix de l’immobilier par EPCI
# On situe le dossier dans lequel se trouve nos données
csv_path <- here("data", "donnees_standr.csv")
# lecture du CSV dans un dataframe
immo_df <- read.csv2(csv_path)
# Pour visualiser les 10 1ères lignes
datatable(head(immo_df, 10))Ce fichier est composé des 10 variables suivantes :
- SIREN : code SIREN de l’EPCI
- prix_med : prix médian par EPCI à la vente (au m²)
- perc_log_vac : % logements vacants
- perc_maison : % maisons
- perc_tiny_log : % petits logements (surface < ?)
- dens_pop : densité de population (nb habitants / km² ?)
- med_niveau_vis : médiane du niveau de vie
- part_log_suroccup : % logements suroccupés
- part_agri_nb_emploi : % agriculteurs
- part_cadre_profintellec_nbemploi : % cadres et professions intellectuelles
La variable SIREN nous servira de “clé” pour joindre ces
données statistiques aux données spatiales, la variable
prix_med sera la variable que nous chercherons à expliquer
(VD), et toutes les autres seront nos variables explicatives (VI).
1.2.2 Chargement des données géographiques : les EPCI de France métropolitaine
Ces données proviennent de la base ADMIN-EPXRESS-COG de l’IGN, édition 2022. Le format d’entrée est le shapefile mais nous passerons par une conversion au format sf, ce qui nous permet d’utiliser le package mapsf, pour les prévisualiser :
# lecture du shapefile en entrée dans un objet sf
shp_path <- here("data", "EPCI.shp")
epci_sf <- st_read(shp_path)
## Reading layer `EPCI' from data source
## `/media/glecampion/Data/gregoire/atelierR_foss_gwr/data/EPCI.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1242 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 124110 ymin: 6049642 xmax: 1242428 ymax: 7110453
## Projected CRS: RGF93 v1 / Lambert-93
# visualisation des données géographiques
mf_map(x = epci_sf)1.2.3 Jointure des données géographiques et tabulaires
Les lignes vides ne nous intéressent pas, nous ne les conserverons pas car elles pourraient poser problème lors de la réalisation de nos analyses. On fait la jointure en ne gardant que les EPCI ayant une correspondance dans le tableau de données :
## [1] 1223
Pourquoi certains EPCI n’ont-ils pas de correspondance dans notre tableau de données ? Il s’agit notamment des EPCI de la petite couronne de Paris, qui se superposent à l’EPCI de la métropole du grand Paris (les données sont donc bien présentes pour cette zone). Pour le reste, il s’agit de nouveaux EPCI ou bien d’EPCI ayant évolué et donc changé d’identifiant.
On peut maintenant représenter les données du prix médian de l’immobilier par EPCI sous forme de carte :
mf_map(x = data_immo,
var = "prix_med",
type = "choro",
breaks = "quantile",
nbreaks = 7,
pal = "Mint",
lwd = 0.01,
leg_title = "Discrétisation par quantile",
leg_val_rnd = 0)
mf_title("Prix médian de l'immobilier au m² par EPCI")
mf_credits("Sources: Notaires de France, IGN Admin Express")Et avoir un aperçu rapide des autres données issues de l’INSEE :
par(mfrow = c(3,3))
for (var in colnames(data_immo)[6:13]) {
mf_map(x = data_immo,
var = var,
type = "choro",
breaks = "quantile",
nbreaks = 7,
pal = "Purples",
lwd = 0.01,
leg_pos = NA)
mf_title(var)
mf_credits('Sources: INSEE, IGN Admin Express')
}Dans le cas où vous préféreriez manipuler vos données sous un format sp (package sp), ou dans le cas où ce format serait requis pour utiliser certains packages ou certaines formules, vous pouvez convertir votre objet sf en sp à l’aide de la ligne de commande suivante (nous en aurons besoin pour la suite) :
Les packages sf (Simple Features for R) et sp (Classes and methods for spatial data) proposent tous deux des fonctions pour manipuler des données spatiales. Le package sf est plus récent et plus performant, mais beaucoup de packages R ne fonctionnent encore qu’avec le format sp.
2 Régression géographiquement pondérée (GWR)
Avant d’attaquer la GWR un petit appel des épisodes précédents !
Notre objectif est donc d’expliquer le prix median de l’immobilier des EPCI en France métropolitaine, en fonction d’un certian nombre d’autre variables comme par exemple la densité de population, le pourcentage de petits logement… Alors que nous pensions pouvoir réaliser une régression classique nous nous sommes rendu compte que cela ne fonctionnait pas ! Le spatial comme à son habitude venait perturbé nos plans.
Avant de poursuivre prenons le temps de visualiter la distribution de nos données le prix de l’immobilier et ses prédicteurs :
# Distribution des variables indépendantes :
a <- add_histogram(plot_ly(data_immo, x = ~log(perc_log_vac), name = "perc_log_vac"))
b <- add_histogram(plot_ly(data_immo, x = ~log(perc_maison), name = "perc_maison"))
c <- add_histogram(plot_ly(data_immo, x = ~log(perc_tiny_log), name = "perc_tiny_log"))
d <- add_histogram(plot_ly(data_immo, x = ~log(dens_pop), name = "dens_pop"))
e <- add_histogram(plot_ly(data_immo, x = ~log(med_niveau_vis), name = "med_niveau_vis"))
f <- add_histogram(plot_ly(data_immo, x = ~log(part_log_suroccup), name = "part_log_suroccup"))
g <- add_histogram(plot_ly(data_immo, x = ~log(part_agri_nb_emploi), name = "part_agri_nb_emploi"))
h <- add_histogram(plot_ly(data_immo, x = ~log(part_cadre_profintellec_nbemploi), name = "part_cadre_profintellec_nbemploi"))
fig = subplot(a, b, c, d, e, f, g, h, nrows = 2)
figVoici à quoi ressemble notre modèle linéaire classique
# Dans le fonctionnement sur R il est important de stocker la régression dans un objet.
# Pour lancer la régression on va utiliser la fonction lm() dont les 2 lettres sont l'acronyme pour linear model
mod.lm <- lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi,
data = data_immo)
# On affiche les principaux résultats avec la fonction summary
summary(mod.lm)##
## Call:
## lm(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log +
## dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi +
## part_cadre_profintellec_nbemploi, data = data_immo)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1566.8 -220.2 -27.2 174.0 3277.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1592.873 11.204 142.171 < 2e-16 ***
## perc_log_vac -287.337 14.134 -20.330 < 2e-16 ***
## perc_maison -128.255 20.041 -6.400 2.22e-10 ***
## perc_tiny_log -261.556 27.934 -9.363 < 2e-16 ***
## dens_pop 173.942 15.272 11.389 < 2e-16 ***
## med_niveau_vis 288.017 13.860 20.781 < 2e-16 ***
## part_log_suroccup 388.982 27.352 14.221 < 2e-16 ***
## part_agri_nb_emploi -20.785 15.059 -1.380 0.168
## part_cadre_profintellec_nbemploi -7.841 19.904 -0.394 0.694
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 391.8 on 1214 degrees of freedom
## Multiple R-squared: 0.7784, Adjusted R-squared: 0.7769
## F-statistic: 532.9 on 8 and 1214 DF, p-value: < 2.2e-16
| Characteristic | Beta | 95% CI1 | p-value |
|---|---|---|---|
| (Intercept) | 1,593 | 1,571, 1,615 | <0.001 |
| perc_log_vac | -287 | -315, -260 | <0.001 |
| perc_maison | -128 | -168, -89 | <0.001 |
| perc_tiny_log | -262 | -316, -207 | <0.001 |
| dens_pop | 174 | 144, 204 | <0.001 |
| med_niveau_vis | 288 | 261, 315 | <0.001 |
| part_log_suroccup | 389 | 335, 443 | <0.001 |
| part_agri_nb_emploi | -21 | -50, 8.8 | 0.2 |
| part_cadre_profintellec_nbemploi | -7.8 | -47, 31 | 0.7 |
| 1 CI = Confidence Interval | |||
On peut également visualiser graphiquement les coefficients des
variables explicatives avec le package GGally :
Une manière très efficace de vérifier l’importance de la dimension spatiale c’est de cartographier nos résidus. S’il semble émerger une structure c’est que très vraissemblablement le spatial va jouer et que les conditions liées aux résidus de la réalisation d’une régression classique ne sont pas vérifiées.
De manière générale lorsque l’on travaille avec des données spatiales il est toujours bon de les cartographier.
On intègre les résidus à la table attributaire de notre objet sf. A priori, comme on a utilisé nos données spatiales (sf) pour la régression les données sont classées dans le bon ordre.
On peut maintenant faire la carte des résidus :
# on fonce légèrement la couleur de fond pour mieux voir la classe centrale
mf_theme("default", bg = "#dedede")
# création d'une palette divergente avec une couleur neutre centrale
palette = mf_get_pal(n = c(nb_cl_inf0_res, nb_cl_sup0_res), pal = c("Teal", "Peach"), neutral = "#f5f5f5")
# la carte :
mf_map(x = data_immo,
var = "res_reg",
type = "choro",
border = NA,
lwd = 0.1,
breaks = breaks_residus,
pal = palette,
leg_title = "Valeur centrale = 0\nIntervalle = σ / 2",
leg_val_rnd = 1)
mf_title("Résidus de régression linéaire classique")
mf_credits("Sources: Notaires de France, INSEE, IGN Admin Express")Sur cette carte on voit très clairement une spatialisation des résidus, sans même faire les tests nous aurions pu voir que la dimension spatiale jouait bien un rôle. Sans autocorrelation nous aurions eu une répartition aléatoire des résidus.
Comme nous l’avons vu, les données, lorsqu’elles sont spatialisées, sont souvent soumises aux phénomènes de dépendance et d’hétérogénéité spatiale.
Pour prendre en compte le phénomène de dépendance spatiale, on fait souvent appel aux régressions spatiales. Elles permettent à la fois de mieux comprendre la relation qui unit les variables explicatives à la variable étudiée d’une part ; et d’autre part de traiter la relation de notre VD à son propre voisinage. Les régressions spatiales peuvent donc (en fonction du type de régression choisie) intégrer les caractéristiques du voisinage (de la VD ou des VI) pour expliquer la VD. Il existe différents modèles de régression spatiale, toute la question est de savoir quel modèle utiliser ? Ce choix va dépendre de la nature des phénomènes spatiaux étudiés. Nous ne traiterons pas ici de ce sujet.
Par ailleurs, l’hétérogénéité, qui renvoie à une instabilité, induit une variabilité spatiale de nos paramètres. L’idée est que nos VI peuvent avoir un effet qui n’est pas le même partout dans l’espace. Dans ce cas nous optons pour la régression géographiquement pondérée (GWR).
Ainsi, si on s’intéresse au lien entre les voisins on est dans l’autocorrélation spatiale et les modèle SEM, SAR & cie mais si on s’intéresse à l’hétérogénéité de nos variables c’est-à-dire à leur variabilité selon leur localisation on est dans la GWR.
Attention, ce qui en théorie peut paraître assez tranché ne l’est souvent pas du tout en pratique. En effet, il y a bien des cas où l’on a du mal à savoir dans quel cadre on se situe exactement.
Pour réaliser une GWR sur R plusieurs packages reconnus existent. On
peut citer notamment le package spgwr
et le package GWmodel.
Nous choisirons d’utiliser ici le package GWmodel.
2.1 Calcul de la matrice des distance
La première étape est de calculer la distance entre toutes nos
observations. Pour ce faire nous utiliserons la fonction
gw.dist().
2.2 Définition de la bande passante
La bande passante est une distance au-delà de laquelle le
poids des observations est considéré comme nul. Le calcul de
cette distance est très important car la valeur de la bande passante
pourra fortement influencer notre modèle. La définition de la bande
passante renvoie à quel type de pondération nous souhaitons appliquer.
Heureusement la fonction bw.gwr va choisir pour nous le
résultat optimal…
Pour ce faire la fonction va se baser sur un critère statistique que l’utilisateur devra définir : le CV (validation croisée) ou le AIC (Critère d’information d’Akaike). Elle reposera aussi sur un type de noyau qu’il faudra également définir : Gaussien, Exponentiel, Bicarré, Tricube ou encore Boxcar. Enfin, il sera également nécessaire de savoir si ce noyau pourra être adaptatif ou fixe.
Voici quelques informations pour guider nos choix :
- Le critère CV a pour objectif de maximiser le pouvoir prédictif du modèle, le critère AIC va chercher un compromis entre le pouvoir prédictif du modèle et son degré de complexité. En général, le critère AIC est privilégié.
- Avec un noyau fixe l’étendue du noyau est déterminée par la distance au point d’intérêt et il est identique en tout point de l’espace. Un noyau fixe est adapté si la répartition des données est homogène dans l’espace, l’unité de la bande passante sera donc une distance. Avec un noyau adaptatif l’étendue du noyau est déterminée par le nombre de voisins. Il est donc plus adapté à une répartition non homogène, l’unité sera alors le nombre de voisins.
Concernant la forme des noyaux :
- Les noyaux gaussiens et exponentiels vont pondérer toutes les observations avec un poids qui tend vers zéro avec la distance au point estimé.
- Les noyaux bisquare et tricube (dont les formes sont très proches) accordent également aux observations un poids décroissant avec la distance, mais par contre ce poids est nul au delà de la distance définie par la bande passante.
- Le noyau Box-Car traite un phénomène continu de façon discontinue.
Manuel de géographie quantitative (Feuillet et al. 2019)
Sachant que sur la forme du noyau, il est tout à fait possible de comparer deux pondérations et deux modèles de GWR.
Ici, nous allons tester avec un noyau gaussien, ce qui sera justifié un petit peu plus bas.
# Définition de la bande passante (bandwidth en anglais) :
bw_g <- bw.gwr(data = data_immo_sp,
approach = "AICc",
kernel = "gaussian",
adaptive = TRUE,
dMat = dm.calib,
formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)## Adaptive bandwidth (number of nearest neighbours): 763 AICc value: 17959.04
## Adaptive bandwidth (number of nearest neighbours): 480 AICc value: 17884.14
## Adaptive bandwidth (number of nearest neighbours): 303 AICc value: 17798.07
## Adaptive bandwidth (number of nearest neighbours): 196 AICc value: 17715.24
## Adaptive bandwidth (number of nearest neighbours): 127 AICc value: 17622.76
## Adaptive bandwidth (number of nearest neighbours): 87 AICc value: 17526.89
## Adaptive bandwidth (number of nearest neighbours): 59 AICc value: 17422.28
## Adaptive bandwidth (number of nearest neighbours): 45 AICc value: 17353.68
## Adaptive bandwidth (number of nearest neighbours): 33 AICc value: 17283.09
## Adaptive bandwidth (number of nearest neighbours): 29 AICc value: 17261.83
## Adaptive bandwidth (number of nearest neighbours): 23 AICc value: 17250.23
## Adaptive bandwidth (number of nearest neighbours): 22 AICc value: 17247.23
## Adaptive bandwidth (number of nearest neighbours): 19 AICc value: 17201.79
## Adaptive bandwidth (number of nearest neighbours): 19 AICc value: 17201.79
## [1] 19
Notre bande passante est donc ici de 19 voisins, ce qui implique que les EPCI au-delà de cette “distance” auront un poids ramené à 0 et ne joueront donc plus de rôle dans la description de la relation statistique.
2.3 Estimation du modèle
Une fois la bande passante définie on peut lancer l’estimation de notre modèle de GWR :
mod.gwr_g <- gwr.robust(data = data_immo_sp,
dMat = dm.calib,
bw = bw_g,
kernel = "gaussian",
filtered = FALSE, # un des problèmes de la GWR est de gérer des individus "aberrants" au niveau local. 2 méthodes ont été définies pour gérer cela.
# Méthode 1 (argument TRUE) on filtre en fonction des individus standardisés. L'objectif est de détecter les individus dont les résidus sont très élevés et de les exclure.
# Méthode 2 (argument FALSE) on diminue le poids des observations aux résidus élevés.
adaptive = TRUE,
formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)Si on souhaite comparer deux modèles car nous avons un doute sur les paramètres c’est tout à fait possible. Par exemple ici nous souhaitons comparer deux formes de noyau :
bw_tri <- bw.gwr(data = data_immo_sp,
approach = "AICc",
kernel = "tricube",
adaptive = TRUE,
dMat = dm.calib,
formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)## Adaptive bandwidth (number of nearest neighbours): 763 AICc value: 17701.17
## Adaptive bandwidth (number of nearest neighbours): 480 AICc value: 17586.94
## Adaptive bandwidth (number of nearest neighbours): 303 AICc value: 17422.51
## Adaptive bandwidth (number of nearest neighbours): 196 AICc value: 17294.75
## Adaptive bandwidth (number of nearest neighbours): 127 AICc value: 17254.1
## Adaptive bandwidth (number of nearest neighbours): 87 AICc value: 17279.33
## Adaptive bandwidth (number of nearest neighbours): 154 AICc value: 17259.05
## Adaptive bandwidth (number of nearest neighbours): 112 AICc value: 17257.17
## Adaptive bandwidth (number of nearest neighbours): 138 AICc value: 17254.8
## Adaptive bandwidth (number of nearest neighbours): 122 AICc value: 17253.75
## Adaptive bandwidth (number of nearest neighbours): 117 AICc value: 17255.04
## Adaptive bandwidth (number of nearest neighbours): 123 AICc value: 17253.57
## Adaptive bandwidth (number of nearest neighbours): 126 AICc value: 17254.25
## Adaptive bandwidth (number of nearest neighbours): 123 AICc value: 17253.57
mod.gwr_tri <- gwr.robust(data = data_immo_sp,
dMat = dm.calib,
bw = bw_tri,
kernel = "gaussian",
filtered = FALSE,
adaptive = TRUE,
formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi)
Best_gwr <- cbind(
rbind(bw_g, bw_tri),
rbind(mod.gwr_g$GW.diagnostic$gw.R2,mod.gwr_tri$GW.diagnostic$gw.R2),
rbind(mod.gwr_g$GW.diagnostic$AIC,mod.gwr_tri$GW.diagnostic$AIC)) %>%
`colnames<-`(c("Nb Voisins","R2","AIC")) %>%
`rownames<-`(c("GAUSSIAN","TRICUBE"))
Best_gwr## Nb Voisins R2 AIC
## GAUSSIAN 19 0.9324328 16849.26
## TRICUBE 123 0.8577370 17567.05
Le modèle avec une forme qui a été définie au format gaussien explique un meilleur \(R^2\) et le score d’\(AIC\) est plus faible. Ce modèle est donc plus performant et dans notre situation c’est plutôt ce modèle qu’il faut privilégier.
2.4 Interprétation des premiers résultats
Comme pour le modèle linéaire classique, l’objet qui contient notre GWR est composé de plusieurs éléments. Pour obtenir nos résultats il suffit d’appeler l’objet.
## Length Class Mode
## GW.arguments 11 -none- list
## GW.diagnostic 8 -none- list
## lm 14 lm list
## SDF 1223 SpatialPolygonsDataFrame S4
## timings 5 -none- list
## this.call 13 -none- call
## Ftests 0 -none- list
## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2023-10-07 14:25:49.656371
## Call:
## gwr.basic(formula = formula, data = data, bw = bw, kernel = kernel,
## adaptive = adaptive, p = p, theta = theta, longlat = longlat,
## dMat = dMat, F123.test = F123.test, cv = T, W.vect = W.vect)
##
## Dependent (y) variable: prix_med
## Independent variables: perc_log_vac perc_maison perc_tiny_log dens_pop med_niveau_vis part_log_suroccup part_agri_nb_emploi part_cadre_profintellec_nbemploi
## Number of data points: 1223
## ***********************************************************************
## * Results of Global Regression *
## ***********************************************************************
##
## Call:
## lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1566.8 -220.2 -27.2 174.0 3277.3
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1592.873 11.204 142.171 < 2e-16 ***
## perc_log_vac -287.337 14.134 -20.330 < 2e-16 ***
## perc_maison -128.255 20.041 -6.400 2.22e-10 ***
## perc_tiny_log -261.556 27.934 -9.363 < 2e-16 ***
## dens_pop 173.942 15.272 11.389 < 2e-16 ***
## med_niveau_vis 288.017 13.860 20.781 < 2e-16 ***
## part_log_suroccup 388.982 27.352 14.221 < 2e-16 ***
## part_agri_nb_emploi -20.785 15.059 -1.380 0.168
## part_cadre_profintellec_nbemploi -7.841 19.904 -0.394 0.694
##
## ---Significance stars
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 391.8 on 1214 degrees of freedom
## Multiple R-squared: 0.7784
## Adjusted R-squared: 0.7769
## F-statistic: 532.9 on 8 and 1214 DF, p-value: < 2.2e-16
## ***Extra Diagnostic information
## Residual sum of squares: 186354789
## Sigma(hat): 390.6721
## AIC: 18086.13
## AICc: 18086.31
## BIC: 16985.31
## ***********************************************************************
## * Results of Geographically Weighted Regression *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: gaussian
## Adaptive bandwidth: 19 (number of nearest neighbours)
## Regression points: the same locations as observations are used.
## Distance metric: A distance matrix is specified for this model calibration.
##
## ****************Summary of GWR coefficient estimates:******************
## Min. 1st Qu. Median
## Intercept 1109.13394 1309.99018 1516.79508
## perc_log_vac -919.82407 -224.60028 -163.40316
## perc_maison -898.16365 -258.64481 -110.20592
## perc_tiny_log -1210.96882 -206.68578 -92.91874
## dens_pop -411.20172 72.82448 217.03593
## med_niveau_vis 81.29015 287.78992 358.32914
## part_log_suroccup -543.68600 42.42839 142.48859
## part_agri_nb_emploi -498.74706 -65.81443 -32.88823
## part_cadre_profintellec_nbemploi -347.37935 -48.57329 0.58811
## 3rd Qu. Max.
## Intercept 1696.07367 2330.78
## perc_log_vac -113.83881 388.64
## perc_maison 76.01079 896.95
## perc_tiny_log 32.81544 773.30
## dens_pop 268.35651 659.84
## med_niveau_vis 413.36560 853.98
## part_log_suroccup 247.33650 700.50
## part_agri_nb_emploi -1.15602 120.33
## part_cadre_profintellec_nbemploi 49.93194 388.50
## ************************Diagnostic information*************************
## Number of data points: 1223
## Effective number of parameters (2trace(S) - trace(S'S)): 320.246
## Effective degrees of freedom (n-2trace(S) + trace(S'S)): 902.754
## AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 17201.79
## AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 16849.26
## BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 17068.01
## Residual sum of squares: 56808914
## R-square value: 0.9324328
## Adjusted R-square value: 0.9084372
##
## ***********************************************************************
## Program stops at: 2023-10-07 14:26:25.464764
Cette visualisation des résultats nous propose d’abord un rappel complet du modèle linéaire classique. Puis viennent ensuite les informations concernant notre GWR. Le premier indicateur à analyser est le \(R^2 ajusté\) ajusté de la GWR, qui est nettement meilleur que celui de la régression linéaire multiple. On passe de \(77\%\) de variance expliqué à \(91\%\) avec la GWR. Ce \(R^2\), pour la GWR, est en fait une moyenne calculée des \(R^2\) de toutes les régressions locales réalisées par la GWR.
La seconde information qui nous intéresse particulièrement est les coefficient associés à nos VI. Nous voyons qu’ils ne sont pas présentés de la même manière que ceux de la régression linéaire. En effet, chaque VI va avoir des coefficients en fonction du minimum, maximum et des quartiles. Cela permet de rendre compte de la stationnarité de l’effet ou non. Dans notre cas on voit qu’il y a bien une variation et même dans certains cas une inversion des signes. Cela laisse supposer une non stationnarité des effets : un effet local peut être présent qui ne suivrait pas l’effet global.
Par exemple, pour le pourcentage de logement vacant avec un coefficient global (modèle linéaire) de \(-287\), quand ce pourcentage augmente le prix médian baisse. En simplifiant le pourcentage baisse d’une unité le prix médian augmente de \(287€\). Dans le cas de la densité de population on a un coefficient global positif donc une relation positive. La densité augmente donc le prix médian augmente. Ici, au global la densité augmente d’une unité le prix médian augmente de \(173€\).
Les résultats de la GWR peuvent donc être lus à une échelle globale pour mesurer la pertinence du modèle ; mais également à des échelles locales : les résultats illustrent ainsi comment les coefficients varient en fonction des unités spatiales. Gardons l’exemple de la densité de population. Dans les lieux où le prix médian est à son minimum le coefficient est de \(-411\) ; on a donc une relation négative. Dans ces espaces la densité augmente d’une unité le prix médian baisse de \(411€\).
Ensuite nous pouvons constater une inversion du signe du coefficient. Ainsi dans les EPCI du dernier quartile où le prix médian du logement est le plus élevé (par ex. Paris) le coefficient est positif. À son maximum une augmentation d’une unité entraîne une augmentation du prix de \(659€\). On a donc très clairement un effet de la densité qui ne sera pas du tout le même en fonction du lieu.
Nous pouvons également étudier l’intervalle interquartile. Ainsi, toujours pour la densité, ce résultat implique que pour 50% de nos unités spatiales (EPCI entre quartile 1 et 3), une augmentation d’une unité de la densité va impliquer une augmentation du prix médian entre \(72€\) et \(268€\).
Au travers de ces résultats on voit parfaitement comment une même variable peut avoir un effet différent, voire opposé en fonction des unités de lieu.
La cartographie va être la meilleure manière de représenter les betas (coefficients) et les différents indicateurs fournis avec la GWR, cela nous permet de décrire plus finement et plus précisément les phénomènes observés.
L’ensemble des données est stocké dans le sous objet SDF de notre modèle. Il contient l’ensemble des informations du modèle associé à chaque donnée spatiale.
On peut le convertir en un dataframe pour le visualiser plus facilement. À l’origine il est au format “SpatialPointsDataFrame”.
# Pour visualiser ce fichier dans R
#View(mod.gwr_g$SDF@data)
#Pour voir à quoi il ressemble dans ce document
datatable(mod.gwr_g$SDF@data)## [1] "Intercept" "perc_log_vac"
## [3] "perc_maison" "perc_tiny_log"
## [5] "dens_pop" "med_niveau_vis"
## [7] "part_log_suroccup" "part_agri_nb_emploi"
## [9] "part_cadre_profintellec_nbemploi" "y"
## [11] "yhat" "residual"
## [13] "CV_Score" "Stud_residual"
## [15] "Intercept_SE" "perc_log_vac_SE"
## [17] "perc_maison_SE" "perc_tiny_log_SE"
## [19] "dens_pop_SE" "med_niveau_vis_SE"
## [21] "part_log_suroccup_SE" "part_agri_nb_emploi_SE"
## [23] "part_cadre_profintellec_nbemploi_SE" "Intercept_TV"
## [25] "perc_log_vac_TV" "perc_maison_TV"
## [27] "perc_tiny_log_TV" "dens_pop_TV"
## [29] "med_niveau_vis_TV" "part_log_suroccup_TV"
## [31] "part_agri_nb_emploi_TV" "part_cadre_profintellec_nbemploi_TV"
## [33] "E_weigts" "Local_R2"
# Intercept : c'est la constante c'est à dire prix médian de référence
# nom de la variable : estimation du coefficient, du beta associé à la VI en chaque point.
# y : les valeurs de la VD
# yhat : valeur de y prédite.
# residual, Stud_residual : résidu et résidu standardisé
# CV_score : score de validation croisée
# _SE : erreur standard de l'estimation du coefficient
# _TV : t-value de l'estimation du coefficient
# E_weight : poids des observations dans la régression robuste
# Local_R2 : R2 au niveau de chaque unité spatiale2.4.1 Étude des résidus
Commençons par une étude des résidus afin de vérifier que cette fois ils n’ont pas de structure apparente.
# on récupère les résidus dans data_immo
res_gwr <- mod.gwr_g$SDF$Stud_residual
data_immo$res_gwr <- res_gwr
# calcul des limites de classes avec la fonction discr, centrées sur 0
res_resgwr <- discr(data_immo$res_gwr, 0, "class_center", sd(data_immo$res_gwr)*0.5, 10)
breaks_gwr <- res_resgwr[[1]]
nb_cl_sup0_gwr <- res_resgwr[[2]]
nb_cl_inf0_gwr <- res_resgwr[[3]]
# création de la palette correspondante
palette = mf_get_pal(n = c(nb_cl_inf0_gwr, nb_cl_sup0_gwr), pal = c("Teal", "Peach"), neutral = "#f5f5f5")
# la carte des résidus
mf_map(x = data_immo,
var = "res_gwr",
type = "choro",
border = "gray",
lwd = 0.2,
breaks = breaks_gwr,
pal = palette,
leg_title = "Discrétisation standardisée :\nvaleur centrale = 0\nintervalle = σ / 2",
leg_val_rnd = 1)
mf_title("Résidus GWR")
mf_credits("Sources: Notaires de France, INSEE, IGN Admin Express")Cette carte ne présente pas de structure spatiale marquée et nous amène à penser que nous avons expliqué l’ensemble des phénomènes spatiaux liés aux questions de prix de l’immobilier.
2.4.2 Étude des coefficients
Pour visualiser la non stationnarité des effets de nos VI la solution la plus efficace est la carte.
# On ajoute à data_immo les coefficients
data_immo$agri.coef <- mod.gwr_g$SDF$part_agri_nb_emploi
data_immo$perc_maison.coef <- mod.gwr_g$SDF$perc_maison
data_immo$dens_pop.coef <- mod.gwr_g$SDF$dens_pop
data_immo$med_vie.coef <- mod.gwr_g$SDF$med_niveau_vis
data_immo$logvac.coef <- mod.gwr_g$SDF$perc_log_vac
data_immo$tinylog.coef <- mod.gwr_g$SDF$perc_tiny_log
data_immo$suroccup.coef <- mod.gwr_g$SDF$part_log_suroccup
data_immo$cadre.coef <- mod.gwr_g$SDF$part_cadre_profintellec_nbemploiPour réaliser les cartes, avec une discrétisation standardisée centrée sur 0 :
par(mfrow = c(4, 2))
for (var in colnames(data_immo)[17:24]) {
# calcul des limites de classe
res <- discr(data.frame(data_immo)[, var], 0, "class_center", sd(data.frame(data_immo)[, var])*0.5, 10)
breaks <- res[[1]]
# palette de couleurs
nb_cl_sup0 <- res[[2]]
nb_cl_inf0 <- res[[3]]
if (nb_cl_inf0 > 0) {
palette = mf_get_pal(n = c(nb_cl_inf0, nb_cl_sup0), pal = c("Teal", "Peach"), neutral = "#f5f5f5")
} else { # cas de la médiane du niveau de vie où la valeur min est supérieure à 0
palette = mf_get_pal(n = c(nb_cl_sup0), pal = c("Peach"))
}
# la carte
mf_map(x = data_immo,
var = var,
type = "choro",
border = "gray",
lwd = 0.1,
breaks = breaks,
pal = palette,
leg_pos = "left",
leg_title = NA,
leg_val_rnd = 0)
mf_title(var)
}Les cartes des betas vont illustrer la variation des effets en fonctions des entités spatiales et de leur voisinage. Dans notre cas on verra quels sont les EPCI où l’effet du coefficient est négatif et ceux où il est positif, c’est-à-dire dans quel EPCI notre VI va entraîner une augmentation du prix médian et dans quel autre au contraire une diminution, toutes choses égales par ailleurs. Sachant que dans notre cas toutes les VI sont significatives, elles ont donc toutes un effet qui varie localement.
Au-delà de cette visualisation VI par VI, il peut être intéressant de savoir par EPCI quelle variable sera la plus explicative dans la relation à notre VD, laquelle a l’impact le plus important. Nous avons donc réalisé une carte des contributions max par EPCI. Pour la réaliser nous nous sommes basés sur le T-value.
data_immo$agri.t <- mod.gwr_g$SDF$part_agri_nb_emploi_TV
data_immo$maison.t <- mod.gwr_g$SDF$perc_maison_TV
data_immo$dens.t <- mod.gwr_g$SDF$dens_pop_TV
data_immo$medvie.t <- mod.gwr_g$SDF$med_niveau_vis_TV
data_immo$logvac.t <- mod.gwr_g$SDF$perc_log_vac_TV
data_immo$tinylog.t <- mod.gwr_g$SDF$perc_tiny_log_TV
data_immo$suroccup.t <- mod.gwr_g$SDF$part_log_suroccup_TV
data_immo$cadre.t <- mod.gwr_g$SDF$part_cadre_profintellec_nbemploi_TV
# Définir contrib max
df <- as.data.frame(data_immo)
# On passe les t-values en valeurs absolues pour voir la plus grande contribution dans un sens sens ou dans l'autre
data_immo$contribmax<- colnames(df[, c(25:32)])[max.col(abs(df[, c(25:32)]),ties.method="first")]par(mfrow = c(1, 1))
# Carte
mf_map(x = data_immo,
var = "contribmax",
type = "typo",
pal = brewer.pal(6,'Set2'),
border = "white",
lwd = 0.2)
mf_title("Carte des variables contribuant le plus par epci")Au-delà de la non-stationnarité de nos VI, cette carte met en évidence un autre phénomène : Dans un modèle linéaire classique, la sélection des VI se fait de manière itérative (ascendante, descendante, ou mixte). L’inclusion d’une VI a pour conséquence d’en exclure d’autres qui présenteraient trop de multi-colinéarité. Le phénomène mis en évidence ici repose sur le fait qu’en fonction du lieu, la hiérarchie des VI, et donc leur pertinence au sein du modèle, n’est pas constante. Cette visualisation des données ouvre donc une perspective intéressante pour la suite. Il serait tout à fait pertinent de développer une méthode permettant d’effectuer un stepwise préalable à toute GWR valable pour chaque lieu individuellement. Un tel modèle traiterait dans sa totalité le problème de non-stationnarité.
Nous pouvons également cartographier les \(R^2\) locaux, ce qui fournit une indication sur les zones où la variabilité est la mieux expliquée.
data_immo$r2local=mod.gwr_g$SDF$Local_R2
mf_map(x = data_immo,
var = "r2local",
type = "choro",
breaks = "quantile",
nbreaks = 11,
pal= "Reds",
border = "gray",
lwd = 0.2,
leg_title = "Discrétisation par quantile")
mf_title("R² locaux")Nous avons choisi d’utiliser une palette de valeurs large pour représenter au mieux les différences de \(R^2\) locaux. Toutefois, cette carte peut paraître trompeuse ; tous les EPCI obtiennent une explication satisfaisante avec un R² local minimum de 0.68.
À partir des t-value on peut aussi étudier la significativité des effets sur le territoire. On peut ainsi calculer et cartographier un indicateur qui représenterait le nombre de VI dont l’effet est significatif sur chaque unité spatiale. Cela donne une bonne idée de la complexité du phénomène sur un espace donné (en effet sur un EPCI on peut avoir toutes les variables significatives, elle jouent donc sur cet espace toutes un rôles) et souligne l’importance d’avoir une carte par coefficient. Cette carte montre également qu’un modèle parfaitement adaptatif gagnerait en sobriété et donnerait un AIC plus satisfaisant en ne sélectionnant localement que les variables réellement significatives dans la relation.
# Pour rappel si on a plus de 200 individus et le t-value > |1.96| on pourra considérer le coefficient comme significatif au seuil de 0.05 (95% chances que ce ne soit pas dû au hasard)
data_immo$nbsignif_t <- rowSums(abs(df[, c(25:32)]) > 1.96)
mf_map(x = data_immo,
var = "nbsignif_t",
type = "typo",
pal = "Reds",
border = "gray",
lwd = 0.2)
mf_title("Nombre de Betas significatifs par EPCI (t-value)")Il se peut que cela soit plus intéressant d’utiliser les p-value, notamment si vous avez moins de 200 individus.
# Les p-value ne sont pas fournis dans le modèle de la GWR, on pourrait les calculer à partir de t-value et de l'erreur standard mais le package GWmodel propose une fonction pour les obtenir
pvalue <- gwr.t.adjust(mod.gwr_g)
# On ajoute les p-value à notre fichier
data_immo$agri.p <- pvalue$SDF$part_agri_nb_emploi_p
data_immo$maison.p <- pvalue$SDF$perc_maison_p
data_immo$dens.p <- pvalue$SDF$dens_pop_p
data_immo$medvie.p <- pvalue$SDF$med_niveau_vis_p
data_immo$logvac.p <- pvalue$SDF$perc_log_vac_p
data_immo$tinylog.p <- pvalue$SDF$perc_tiny_log_p
data_immo$suroccup.p <- pvalue$SDF$part_log_suroccup_p
data_immo$cadre.p <- pvalue$SDF$part_cadre_profintellec_nbemploi_p
df<- as.data.frame(data_immo)
data_immo$nbsignif_p <- rowSums(df[, c(36:43)] < 0.05)
mf_map(x = data_immo,
var = "nbsignif_p",
type = "typo",
pal= "Reds",
border = "gray",
lwd = 0.2,)
mf_title("Nombre des Betas significatifs par EPCI (p-value)")Sur ces deux dernières cartes (nombre de betas significatifs par EPCI avec les t-value et p-value) nous avons pris volontairement une liberté avec les règles de sémiologie graphique de Bertin. En effet, s’agissant de valeurs discrètes, nous aurions dû les représenter sous forme de symboles proportionnels. Cependant, dans la mesure où les EPCI ont une taille relativement homogène et par souci de visibilité, et aussi parce qu’il s’agit plus d’une carte exploratoire que d’un véritable rendu, nous avons opté pour une carte choroplèthe. Que les mânes de Bertin nous pardonnent !
Dans ce cadre, il est possible de réaliser une collection de cartes des p-value (ou t-value) comme ce qui a été fait pour les coefficients. L’intérêt est de voir où l’effet de la VI est significatif et où il ne l’est pas.
# Ici nous représenterons les p-value avec un découpage par classe de significativité et seulement les p-value de 2 VI
par(mfrow = c(1, 2))
# Par exemple les p-value des coefficients de la variable part de l'emploi agriculteur
data_immo<- data_immo %>% mutate(agri.p_fac = case_when(agri.p<= 0.002 ~ "[0;0.002[",
agri.p <= 0.01 ~ "[0.002;0.01[",
agri.p <= 0.05 ~ "[0.01;0.05[",
agri.p <= 0.1 ~ "[0.05;0.1[",
TRUE ~ "[0.1;1]")) %>%
mutate(agri.p_fac = factor(agri.p_fac,
levels = c("[0;0.002[", "[0.002;0.01[",
"[0.01;0.05[",
"[0.05;0.1[", "[0.1;1]")))
mypal2 <- mf_get_pal(n = 5, palette = "OrRd")
mf_map(x = data_immo,
var = "agri.p_fac",
type = "typo",
border = "grey3",
lwd = 0.1,
pal = mypal2,
leg_title = "Classe P-value")
mf_title("P-value du coefficient de la part d'emploi agriculteurs")
# Pour la densité de population
data_immo<- data_immo %>% mutate(dens.p_fac = case_when(dens.p <= 0.002 ~ "[0;0.002[",
dens.p <= 0.01 ~ "[0.002;0.01[",
dens.p <= 0.05 ~ "[0.01;0.05[",
dens.p <= 0.1 ~ "[0.05;0.1[",
TRUE ~ "[0.1;1]")) %>%
mutate(dens.p_fac = factor(dens.p_fac,
levels = c("[0;0.002[", "[0.002;0.01[",
"[0.01;0.05[",
"[0.05;0.1[", "[0.1;1]")))
mypal2 <- mf_get_pal(n = 5, palette = "OrRd")
mf_map(x = data_immo,
var = "dens.p_fac",
type = "typo",
border = "grey3",
lwd = 0.1,
pal=mypal2,
leg_title = "Classe P-value")
mf_title("P-value du coefficient de la densité de population")Ces cartes des p-value sont particulièrement importantes car elles nous donnent les endroits où l’effet est significatif. En effet, on sait que la VI a effet global qui est significatif, qu’elle a en plus une variabilité locale or localement elle n’est pas partout significative. Pour la part d’agriculteur dans l’emploi, l’effet est significatif quasiment uniquement dans le Sud-Est.
3 Pour aller plus loin…
Cette partie et le code qui va avec est directement tiré du cours donnée à SIGR par Thierry Feuillet (feuillet?).
3.1 GWR Multiscalaire
Selon Thierry Feuillet concernant la GWR multiscalaire :
“Il n’y a pas de raison de penser que tous les prédicteurs agissent sur le prix à la même échelle (c’est-à-dire selon un même schéma de voisinage). Certains processus peuvent être locaux, d’autres globaux. Récemment, une extension de la GWR a été proposée, permettant de relâcher cette hypothèse d’égalité des échelles : la GWR multiscalaire (MGWR, Fotheringham et al., 2017). Le principe est simple : un algorithme optimise le choix de la bandwidth pour chaque prédicteur, en fonction des autres. Il en résulte un modèle souvent mixte.”
Avec notre modèle même simplifié avec 3 prédicteurs le modèle est très long à tourner à l’échelle de la France entière, plus de 2h !
On part donc ici sur un corpus plus petit : La Bretagne sans toucher à nos VI et au nombre d’itérationsLa GWR multiscalaire est une méthode itérative, la faire tourner peut donc être TRES long. Plusieurs possibilités : diminuer le nombre maximum d’itérations, simplifier le modèle et intégrer moins de VI ou sinon réduire l’emprise spatiale.
Afin de réduire les temps de traitement, on filtre d’abord les données pour ne garder que les EPCI bretons. Un EPCI pouvant être à cheval sur 2 régions, on ne va garder ici que les EPCI dont le centroïde est en région Bretagne. On va pour ça se servir de la couche REGION.shp de la base ADMIN-EXPRESS de l’IGN.
# Lecture de la couche région dans un objet sf
shp_path <- here("data", "REGION.shp")
region_sf <- st_read(shp_path)## Reading layer `REGION' from data source
## `/media/glecampion/Data/gregoire/atelierR_foss_gwr/data/REGION.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 13 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 99040 ymin: 6046546 xmax: 1242445 ymax: 7110479
## Projected CRS: RGF93 v1 / Lambert-93
# Pour ne garder que la Bretagne
bzh_sf <- region_sf[region_sf$NOM == "Bretagne",]
mf_map(x = bzh_sf)## Warning: st_centroid assumes attributes are constant over geometries
# Sélection des centroïdes dans la région Bretagne
epci_centroids_bzh <- st_intersection(epci_centroids, bzh_sf)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
# On ne garde que les EPCI correspondant à ces centroïdes
epci_bzh <- merge(x = epci_sf, y = st_drop_geometry(epci_centroids_bzh), by.x = "CODE_SIREN", by.y = "CODE_SIREN")
mf_map(x = epci_bzh)# Il faut refaire la jointure entre les données sur les EPCI et l'objet sf
data_immo_bzh <- merge(x = epci_bzh, y = immo_df, by.x = "CODE_SIREN", by.y = "SIREN")
# conversion objet sf vers objet sp pour le package GWmodel
data_immo_bzh_sp <-as(data_immo_bzh, "Spatial")On peut maintenant lancer la GWR multiscalaire :
#source("gwr.multiscale_T.r")
# On lance la MGWR
# On note qu'il n'est pas ici nécessaire de définir une bande passante. Le principe de la GWR multiscalaire est justement d'adapter à chaque relation locale une bande passante (d'où les itérations)
MGWR <- gwr.multiscale(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi,
data = data_immo_bzh_sp,
kernel = "gaussian",
predictor.centered = rep(T, 3), # centrage des prédicteurs
adaptive = TRUE,
bws0 = rep(1,4)) # BW minimum pour l'optimisation## Warning in gwr.multiscale(formula = prix_med ~ perc_log_vac + perc_maison + :
## All the predictors will be centered, please check the parameter
## predictor.centered
## ------ Calculate the initial beta0 from the above bandwidths ------
## ------ The end for calculating the initial beta0 ------
## ------ Select the optimum bandwidths for each independent variable via AIC aproach ------
## ***** The back-fitting process for model calibration with bandwiths selected *****
## Iteration 1 :
## Now select an optimum bandwidth for the variable: Intercept
## The newly selected bandwidth for variable Intercept is: 14
## The bandwidth used in the last ieration is: 1 and the difference between these two bandwidths is: 13
## The bandwidth for variable Intercept will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: perc_log_vac
## The newly selected bandwidth for variable perc_log_vac is: 53
## The bandwidth used in the last ieration is: 1 and the difference between these two bandwidths is: 52
## The bandwidth for variable perc_log_vac will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: perc_maison
## The newly selected bandwidth for variable perc_maison is: 39
## The bandwidth used in the last ieration is: 1 and the difference between these two bandwidths is: 38
## The bandwidth for variable perc_maison will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: perc_tiny_log
## The newly selected bandwidth for variable perc_tiny_log is: 39
## The bandwidth used in the last ieration is: 1 and the difference between these two bandwidths is: 38
## The bandwidth for variable perc_tiny_log will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: dens_pop
## The newly selected bandwidth for variable dens_pop is: 54
## The bandwidth used in the last ieration is: 1 and the difference between these two bandwidths is: 53
## The bandwidth for variable dens_pop will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: med_niveau_vis
## The newly selected bandwidth for variable med_niveau_vis is: 56
## The bandwidth used in the last ieration is: 1 and the difference between these two bandwidths is: 55
## The bandwidth for variable med_niveau_vis will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: part_log_suroccup
## The newly selected bandwidth for variable part_log_suroccup is: 56
## The bandwidth used in the last ieration is: 1 and the difference between these two bandwidths is: 55
## The bandwidth for variable part_log_suroccup will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: part_agri_nb_emploi
## The newly selected bandwidth for variable part_agri_nb_emploi is: 56
## The bandwidth used in the last ieration is: 1 and the difference between these two bandwidths is: 55
## The bandwidth for variable part_agri_nb_emploi will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: part_cadre_profintellec_nbemploi
## The newly selected bandwidth for variable part_cadre_profintellec_nbemploi is: 56
## The bandwidth used in the last ieration is: 1 and the difference between these two bandwidths is: 55
## The bandwidth for variable part_cadre_profintellec_nbemploi will be continually selected in the next ieration.
## Ieration 1 the differential change value of RSS (dCVR) is: 0.4470374
## ----------End of Iteration 1 ----------
## Iteration 2 :
## Now select an optimum bandwidth for the variable: Intercept
## The newly selected bandwidth for variable Intercept is: 11
## The bandwidth used in the last ieration is: 14 and the difference between these two bandwidths is: 3
## The bandwidth for variable Intercept will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: perc_log_vac
## The newly selected bandwidth for variable perc_log_vac is: 53
## The bandwidth used in the last ieration is: 53 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_log_vac seems to be converged for 1 times.It will be continually optimized in the next 4 times
## Now select an optimum bandwidth for the variable: perc_maison
## The newly selected bandwidth for variable perc_maison is: 56
## The bandwidth used in the last ieration is: 39 and the difference between these two bandwidths is: 17
## The bandwidth for variable perc_maison will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: perc_tiny_log
## The newly selected bandwidth for variable perc_tiny_log is: 39
## The bandwidth used in the last ieration is: 39 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_tiny_log seems to be converged for 1 times.It will be continually optimized in the next 4 times
## Now select an optimum bandwidth for the variable: dens_pop
## The newly selected bandwidth for variable dens_pop is: 56
## The bandwidth used in the last ieration is: 54 and the difference between these two bandwidths is: 2
## The bandwidth for variable dens_pop will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: med_niveau_vis
## The newly selected bandwidth for variable med_niveau_vis is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable med_niveau_vis seems to be converged for 1 times.It will be continually optimized in the next 4 times
## Now select an optimum bandwidth for the variable: part_log_suroccup
## The newly selected bandwidth for variable part_log_suroccup is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_log_suroccup seems to be converged for 1 times.It will be continually optimized in the next 4 times
## Now select an optimum bandwidth for the variable: part_agri_nb_emploi
## The newly selected bandwidth for variable part_agri_nb_emploi is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_agri_nb_emploi seems to be converged for 1 times.It will be continually optimized in the next 4 times
## Now select an optimum bandwidth for the variable: part_cadre_profintellec_nbemploi
## The newly selected bandwidth for variable part_cadre_profintellec_nbemploi is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_cadre_profintellec_nbemploi seems to be converged for 1 times.It will be continually optimized in the next 4 times
## Ieration 2 the differential change value of RSS (dCVR) is: 0.2210537
## ----------End of Iteration 2 ----------
## Iteration 3 :
## Now select an optimum bandwidth for the variable: Intercept
## The newly selected bandwidth for variable Intercept is: 11
## The bandwidth used in the last ieration is: 11 and the difference between these two bandwidths is: 0
## The bandwidth for variable Intercept seems to be converged for 1 times.It will be continually optimized in the next 4 times
## Now select an optimum bandwidth for the variable: perc_log_vac
## The newly selected bandwidth for variable perc_log_vac is: 53
## The bandwidth used in the last ieration is: 53 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_log_vac seems to be converged for 2 times.It will be continually optimized in the next 3 times
## Now select an optimum bandwidth for the variable: perc_maison
## The newly selected bandwidth for variable perc_maison is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_maison seems to be converged for 1 times.It will be continually optimized in the next 4 times
## Now select an optimum bandwidth for the variable: perc_tiny_log
## The newly selected bandwidth for variable perc_tiny_log is: 39
## The bandwidth used in the last ieration is: 39 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_tiny_log seems to be converged for 2 times.It will be continually optimized in the next 3 times
## Now select an optimum bandwidth for the variable: dens_pop
## The newly selected bandwidth for variable dens_pop is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable dens_pop seems to be converged for 1 times.It will be continually optimized in the next 4 times
## Now select an optimum bandwidth for the variable: med_niveau_vis
## The newly selected bandwidth for variable med_niveau_vis is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable med_niveau_vis seems to be converged for 2 times.It will be continually optimized in the next 3 times
## Now select an optimum bandwidth for the variable: part_log_suroccup
## The newly selected bandwidth for variable part_log_suroccup is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_log_suroccup seems to be converged for 2 times.It will be continually optimized in the next 3 times
## Now select an optimum bandwidth for the variable: part_agri_nb_emploi
## The newly selected bandwidth for variable part_agri_nb_emploi is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_agri_nb_emploi seems to be converged for 2 times.It will be continually optimized in the next 3 times
## Now select an optimum bandwidth for the variable: part_cadre_profintellec_nbemploi
## The newly selected bandwidth for variable part_cadre_profintellec_nbemploi is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_cadre_profintellec_nbemploi seems to be converged for 2 times.It will be continually optimized in the next 3 times
## Ieration 3 the differential change value of RSS (dCVR) is: 0.1316818
## ----------End of Iteration 3 ----------
## Iteration 4 :
## Now select an optimum bandwidth for the variable: Intercept
## The newly selected bandwidth for variable Intercept is: 11
## The bandwidth used in the last ieration is: 11 and the difference between these two bandwidths is: 0
## The bandwidth for variable Intercept seems to be converged for 2 times.It will be continually optimized in the next 3 times
## Now select an optimum bandwidth for the variable: perc_log_vac
## The newly selected bandwidth for variable perc_log_vac is: 53
## The bandwidth used in the last ieration is: 53 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_log_vac seems to be converged for 3 times.It will be continually optimized in the next 2 times
## Now select an optimum bandwidth for the variable: perc_maison
## The newly selected bandwidth for variable perc_maison is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_maison seems to be converged for 2 times.It will be continually optimized in the next 3 times
## Now select an optimum bandwidth for the variable: perc_tiny_log
## The newly selected bandwidth for variable perc_tiny_log is: 39
## The bandwidth used in the last ieration is: 39 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_tiny_log seems to be converged for 3 times.It will be continually optimized in the next 2 times
## Now select an optimum bandwidth for the variable: dens_pop
## The newly selected bandwidth for variable dens_pop is: 54
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 2
## The bandwidth for variable dens_pop will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: med_niveau_vis
## The newly selected bandwidth for variable med_niveau_vis is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable med_niveau_vis seems to be converged for 3 times.It will be continually optimized in the next 2 times
## Now select an optimum bandwidth for the variable: part_log_suroccup
## The newly selected bandwidth for variable part_log_suroccup is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_log_suroccup seems to be converged for 3 times.It will be continually optimized in the next 2 times
## Now select an optimum bandwidth for the variable: part_agri_nb_emploi
## The newly selected bandwidth for variable part_agri_nb_emploi is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_agri_nb_emploi seems to be converged for 3 times.It will be continually optimized in the next 2 times
## Now select an optimum bandwidth for the variable: part_cadre_profintellec_nbemploi
## The newly selected bandwidth for variable part_cadre_profintellec_nbemploi is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_cadre_profintellec_nbemploi seems to be converged for 3 times.It will be continually optimized in the next 2 times
## Ieration 4 the differential change value of RSS (dCVR) is: 0.1174028
## ----------End of Iteration 4 ----------
## Iteration 5 :
## Now select an optimum bandwidth for the variable: Intercept
## The newly selected bandwidth for variable Intercept is: 11
## The bandwidth used in the last ieration is: 11 and the difference between these two bandwidths is: 0
## The bandwidth for variable Intercept seems to be converged for 3 times.It will be continually optimized in the next 2 times
## Now select an optimum bandwidth for the variable: perc_log_vac
## The newly selected bandwidth for variable perc_log_vac is: 53
## The bandwidth used in the last ieration is: 53 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_log_vac seems to be converged for 4 times.It will be continually optimized in the next 1 times
## Now select an optimum bandwidth for the variable: perc_maison
## The newly selected bandwidth for variable perc_maison is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_maison seems to be converged for 3 times.It will be continually optimized in the next 2 times
## Now select an optimum bandwidth for the variable: perc_tiny_log
## The newly selected bandwidth for variable perc_tiny_log is: 39
## The bandwidth used in the last ieration is: 39 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_tiny_log seems to be converged for 4 times.It will be continually optimized in the next 1 times
## Now select an optimum bandwidth for the variable: dens_pop
## The newly selected bandwidth for variable dens_pop is: 47
## The bandwidth used in the last ieration is: 54 and the difference between these two bandwidths is: 7
## The bandwidth for variable dens_pop will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: med_niveau_vis
## The newly selected bandwidth for variable med_niveau_vis is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable med_niveau_vis seems to be converged for 4 times.It will be continually optimized in the next 1 times
## Now select an optimum bandwidth for the variable: part_log_suroccup
## The newly selected bandwidth for variable part_log_suroccup is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_log_suroccup seems to be converged for 4 times.It will be continually optimized in the next 1 times
## Now select an optimum bandwidth for the variable: part_agri_nb_emploi
## The newly selected bandwidth for variable part_agri_nb_emploi is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_agri_nb_emploi seems to be converged for 4 times.It will be continually optimized in the next 1 times
## Now select an optimum bandwidth for the variable: part_cadre_profintellec_nbemploi
## The newly selected bandwidth for variable part_cadre_profintellec_nbemploi is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_cadre_profintellec_nbemploi seems to be converged for 4 times.It will be continually optimized in the next 1 times
## Ieration 5 the differential change value of RSS (dCVR) is: 0.1086761
## ----------End of Iteration 5 ----------
## Iteration 6 :
## Now select an optimum bandwidth for the variable: Intercept
## The newly selected bandwidth for variable Intercept is: 11
## The bandwidth used in the last ieration is: 11 and the difference between these two bandwidths is: 0
## The bandwidth for variable Intercept seems to be converged for 4 times.It will be continually optimized in the next 1 times
## Now select an optimum bandwidth for the variable: perc_log_vac
## The newly selected bandwidth for variable perc_log_vac is: 53
## The bandwidth used in the last ieration is: 53 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_log_vac seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable: perc_maison
## The newly selected bandwidth for variable perc_maison is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_maison seems to be converged for 4 times.It will be continually optimized in the next 1 times
## Now select an optimum bandwidth for the variable: perc_tiny_log
## The newly selected bandwidth for variable perc_tiny_log is: 39
## The bandwidth used in the last ieration is: 39 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_tiny_log seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable: dens_pop
## The newly selected bandwidth for variable dens_pop is: 41
## The bandwidth used in the last ieration is: 47 and the difference between these two bandwidths is: 6
## The bandwidth for variable dens_pop will be continually selected in the next ieration.
## Now select an optimum bandwidth for the variable: med_niveau_vis
## The newly selected bandwidth for variable med_niveau_vis is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable med_niveau_vis seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable: part_log_suroccup
## The newly selected bandwidth for variable part_log_suroccup is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_log_suroccup seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable: part_agri_nb_emploi
## The newly selected bandwidth for variable part_agri_nb_emploi is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_agri_nb_emploi seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable: part_cadre_profintellec_nbemploi
## The newly selected bandwidth for variable part_cadre_profintellec_nbemploi is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable part_cadre_profintellec_nbemploi seems to be converged and will be kept the same in the following ierations.
## Ieration 6 the differential change value of RSS (dCVR) is: 0.101719
## ----------End of Iteration 6 ----------
## Iteration 7 :
## Now select an optimum bandwidth for the variable: Intercept
## The newly selected bandwidth for variable Intercept is: 11
## The bandwidth used in the last ieration is: 11 and the difference between these two bandwidths is: 0
## The bandwidth for variable Intercept seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable: perc_maison
## The newly selected bandwidth for variable perc_maison is: 56
## The bandwidth used in the last ieration is: 56 and the difference between these two bandwidths is: 0
## The bandwidth for variable perc_maison seems to be converged and will be kept the same in the following ierations.
## Now select an optimum bandwidth for the variable: dens_pop
## The newly selected bandwidth for variable dens_pop is: 41
## The bandwidth used in the last ieration is: 41 and the difference between these two bandwidths is: 0
## The bandwidth for variable dens_pop seems to be converged for 1 times.It will be continually optimized in the next 4 times
## Ieration 7 the differential change value of RSS (dCVR) is: 0.07537341
## ----------End of Iteration 7 ----------
## Iteration 8 :
## Now select an optimum bandwidth for the variable: dens_pop
## The newly selected bandwidth for variable dens_pop is: 41
## The bandwidth used in the last ieration is: 41 and the difference between these two bandwidths is: 0
## The bandwidth for variable dens_pop seems to be converged for 2 times.It will be continually optimized in the next 3 times
## Ieration 8 the differential change value of RSS (dCVR) is: 0.06691285
## ----------End of Iteration 8 ----------
## Iteration 9 :
## Now select an optimum bandwidth for the variable: dens_pop
## The newly selected bandwidth for variable dens_pop is: 41
## The bandwidth used in the last ieration is: 41 and the difference between these two bandwidths is: 0
## The bandwidth for variable dens_pop seems to be converged for 3 times.It will be continually optimized in the next 2 times
## Ieration 9 the differential change value of RSS (dCVR) is: 0.05956372
## ----------End of Iteration 9 ----------
## Iteration 10 :
## Now select an optimum bandwidth for the variable: dens_pop
## The newly selected bandwidth for variable dens_pop is: 41
## The bandwidth used in the last ieration is: 41 and the difference between these two bandwidths is: 0
## The bandwidth for variable dens_pop seems to be converged for 4 times.It will be continually optimized in the next 1 times
## Ieration 10 the differential change value of RSS (dCVR) is: 0.05253431
## ----------End of Iteration 10 ----------
## Iteration 11 :
## Now select an optimum bandwidth for the variable: dens_pop
## The newly selected bandwidth for variable dens_pop is: 41
## The bandwidth used in the last ieration is: 41 and the difference between these two bandwidths is: 0
## The bandwidth for variable dens_pop seems to be converged and will be kept the same in the following ierations.
## Ieration 11 the differential change value of RSS (dCVR) is: 0.04631346
## ----------End of Iteration 11 ----------
## Iteration 12 :
## Ieration 12 the differential change value of RSS (dCVR) is: 0.04089679
## ----------End of Iteration 12 ----------
## Iteration 13 :
## Ieration 13 the differential change value of RSS (dCVR) is: 0.0361525
## ----------End of Iteration 13 ----------
## Iteration 14 :
## Ieration 14 the differential change value of RSS (dCVR) is: 0.03197006
## ----------End of Iteration 14 ----------
## Iteration 15 :
## Ieration 15 the differential change value of RSS (dCVR) is: 0.02826441
## ----------End of Iteration 15 ----------
## Iteration 16 :
## Ieration 16 the differential change value of RSS (dCVR) is: 0.02496654
## ----------End of Iteration 16 ----------
## Iteration 17 :
## Ieration 17 the differential change value of RSS (dCVR) is: 0.02201899
## ----------End of Iteration 17 ----------
## Iteration 18 :
## Ieration 18 the differential change value of RSS (dCVR) is: 0.01937346
## ----------End of Iteration 18 ----------
## Iteration 19 :
## Ieration 19 the differential change value of RSS (dCVR) is: 0.01698891
## ----------End of Iteration 19 ----------
## Iteration 20 :
## Ieration 20 the differential change value of RSS (dCVR) is: 0.01483008
## ----------End of Iteration 20 ----------
## Iteration 21 :
## Ieration 21 the differential change value of RSS (dCVR) is: 0.01286621
## ----------End of Iteration 21 ----------
## Iteration 22 :
## Ieration 22 the differential change value of RSS (dCVR) is: 0.01106978
## ----------End of Iteration 22 ----------
## Iteration 23 :
## Ieration 23 the differential change value of RSS (dCVR) is: 0.009415249
## ----------End of Iteration 23 ----------
## Iteration 24 :
## Ieration 24 the differential change value of RSS (dCVR) is: 0.007877229
## ----------End of Iteration 24 ----------
## Iteration 25 :
## Ieration 25 the differential change value of RSS (dCVR) is: 0.006427512
## ----------End of Iteration 25 ----------
## Iteration 26 :
## Ieration 26 the differential change value of RSS (dCVR) is: 0.005028265
## ----------End of Iteration 26 ----------
## Iteration 27 :
## Ieration 27 the differential change value of RSS (dCVR) is: 0.003610907
## ----------End of Iteration 27 ----------
## Iteration 28 :
## Ieration 28 the differential change value of RSS (dCVR) is: 0.001957864
## ----------End of Iteration 28 ----------
## Iteration 29 :
## Ieration 29 the differential change value of RSS (dCVR) is: 0.001713618
## ----------End of Iteration 29 ----------
## Iteration 30 :
## Ieration 30 the differential change value of RSS (dCVR) is: 0.00278589
## ----------End of Iteration 30 ----------
## Iteration 31 :
## Ieration 31 the differential change value of RSS (dCVR) is: 0.003322762
## ----------End of Iteration 31 ----------
## Iteration 32 :
## Ieration 32 the differential change value of RSS (dCVR) is: 0.003619623
## ----------End of Iteration 32 ----------
## Iteration 33 :
## Ieration 33 the differential change value of RSS (dCVR) is: 0.003769717
## ----------End of Iteration 33 ----------
## Iteration 34 :
## Ieration 34 the differential change value of RSS (dCVR) is: 0.003819166
## ----------End of Iteration 34 ----------
## Iteration 35 :
## Ieration 35 the differential change value of RSS (dCVR) is: 0.00379544
## ----------End of Iteration 35 ----------
## Iteration 36 :
## Ieration 36 the differential change value of RSS (dCVR) is: 0.003716656
## ----------End of Iteration 36 ----------
## Iteration 37 :
## Ieration 37 the differential change value of RSS (dCVR) is: 0.003595539
## ----------End of Iteration 37 ----------
## Iteration 38 :
## Ieration 38 the differential change value of RSS (dCVR) is: 0.00344137
## ----------End of Iteration 38 ----------
## Iteration 39 :
## Ieration 39 the differential change value of RSS (dCVR) is: 0.003261068
## ----------End of Iteration 39 ----------
## Iteration 40 :
## Ieration 40 the differential change value of RSS (dCVR) is: 0.003059802
## ----------End of Iteration 40 ----------
## Iteration 41 :
## Ieration 41 the differential change value of RSS (dCVR) is: 0.002841355
## ----------End of Iteration 41 ----------
## Iteration 42 :
## Ieration 42 the differential change value of RSS (dCVR) is: 0.002608304
## ----------End of Iteration 42 ----------
## Iteration 43 :
## Ieration 43 the differential change value of RSS (dCVR) is: 0.002362035
## ----------End of Iteration 43 ----------
## Iteration 44 :
## Ieration 44 the differential change value of RSS (dCVR) is: 0.002102565
## ----------End of Iteration 44 ----------
## Iteration 45 :
## Ieration 45 the differential change value of RSS (dCVR) is: 0.001827987
## ----------End of Iteration 45 ----------
## Iteration 46 :
## Ieration 46 the differential change value of RSS (dCVR) is: 0.001533024
## ----------End of Iteration 46 ----------
## Iteration 47 :
## Ieration 47 the differential change value of RSS (dCVR) is: 0.001204747
## ----------End of Iteration 47 ----------
## Iteration 48 :
## Ieration 48 the differential change value of RSS (dCVR) is: 0.0008043307
## ----------End of Iteration 48 ----------
## Iteration 49 :
## Ieration 49 the differential change value of RSS (dCVR) is: 0.000253014
## ----------End of Iteration 49 ----------
## Iteration 50 :
## Ieration 50 the differential change value of RSS (dCVR) is: 0.0008271053
## ----------End of Iteration 50 ----------
## Iteration 51 :
## Ieration 51 the differential change value of RSS (dCVR) is: 0.001103321
## ----------End of Iteration 51 ----------
## Iteration 52 :
## Ieration 52 the differential change value of RSS (dCVR) is: 0.001291758
## ----------End of Iteration 52 ----------
## Iteration 53 :
## Ieration 53 the differential change value of RSS (dCVR) is: 0.001429585
## ----------End of Iteration 53 ----------
## Iteration 54 :
## Ieration 54 the differential change value of RSS (dCVR) is: 0.001532507
## ----------End of Iteration 54 ----------
## Iteration 55 :
## Ieration 55 the differential change value of RSS (dCVR) is: 0.001609177
## ----------End of Iteration 55 ----------
## Iteration 56 :
## Ieration 56 the differential change value of RSS (dCVR) is: 0.001665106
## ----------End of Iteration 56 ----------
## Iteration 57 :
## Ieration 57 the differential change value of RSS (dCVR) is: 0.001704137
## ----------End of Iteration 57 ----------
## Iteration 58 :
## Ieration 58 the differential change value of RSS (dCVR) is: 0.001729128
## ----------End of Iteration 58 ----------
## Iteration 59 :
## Ieration 59 the differential change value of RSS (dCVR) is: 0.001742305
## ----------End of Iteration 59 ----------
## Iteration 60 :
## Ieration 60 the differential change value of RSS (dCVR) is: 0.001745454
## ----------End of Iteration 60 ----------
## Iteration 61 :
## Ieration 61 the differential change value of RSS (dCVR) is: 0.001740053
## ----------End of Iteration 61 ----------
## Iteration 62 :
## Ieration 62 the differential change value of RSS (dCVR) is: 0.001727341
## ----------End of Iteration 62 ----------
## Iteration 63 :
## Ieration 63 the differential change value of RSS (dCVR) is: 0.001708376
## ----------End of Iteration 63 ----------
## Iteration 64 :
## Ieration 64 the differential change value of RSS (dCVR) is: 0.00168407
## ----------End of Iteration 64 ----------
## Iteration 65 :
## Ieration 65 the differential change value of RSS (dCVR) is: 0.001655218
## ----------End of Iteration 65 ----------
## Iteration 66 :
## Ieration 66 the differential change value of RSS (dCVR) is: 0.001622514
## ----------End of Iteration 66 ----------
## Iteration 67 :
## Ieration 67 the differential change value of RSS (dCVR) is: 0.001586571
## ----------End of Iteration 67 ----------
## Iteration 68 :
## Ieration 68 the differential change value of RSS (dCVR) is: 0.001547928
## ----------End of Iteration 68 ----------
## Iteration 69 :
## Ieration 69 the differential change value of RSS (dCVR) is: 0.001507062
## ----------End of Iteration 69 ----------
## Iteration 70 :
## Ieration 70 the differential change value of RSS (dCVR) is: 0.0014644
## ----------End of Iteration 70 ----------
## Iteration 71 :
## Ieration 71 the differential change value of RSS (dCVR) is: 0.001420316
## ----------End of Iteration 71 ----------
## Iteration 72 :
## Ieration 72 the differential change value of RSS (dCVR) is: 0.001375146
## ----------End of Iteration 72 ----------
## Iteration 73 :
## Ieration 73 the differential change value of RSS (dCVR) is: 0.001329187
## ----------End of Iteration 73 ----------
## Iteration 74 :
## Ieration 74 the differential change value of RSS (dCVR) is: 0.001282702
## ----------End of Iteration 74 ----------
## Iteration 75 :
## Ieration 75 the differential change value of RSS (dCVR) is: 0.001235925
## ----------End of Iteration 75 ----------
## Iteration 76 :
## Ieration 76 the differential change value of RSS (dCVR) is: 0.001189062
## ----------End of Iteration 76 ----------
## Iteration 77 :
## Ieration 77 the differential change value of RSS (dCVR) is: 0.001142296
## ----------End of Iteration 77 ----------
## Iteration 78 :
## Ieration 78 the differential change value of RSS (dCVR) is: 0.001095787
## ----------End of Iteration 78 ----------
## Iteration 79 :
## Ieration 79 the differential change value of RSS (dCVR) is: 0.001049675
## ----------End of Iteration 79 ----------
## Iteration 80 :
## Ieration 80 the differential change value of RSS (dCVR) is: 0.001004084
## ----------End of Iteration 80 ----------
## Iteration 81 :
## Ieration 81 the differential change value of RSS (dCVR) is: 0.0009591194
## ----------End of Iteration 81 ----------
## Iteration 82 :
## Ieration 82 the differential change value of RSS (dCVR) is: 0.000914873
## ----------End of Iteration 82 ----------
## Iteration 83 :
## Ieration 83 the differential change value of RSS (dCVR) is: 0.0008714228
## ----------End of Iteration 83 ----------
## Iteration 84 :
## Ieration 84 the differential change value of RSS (dCVR) is: 0.0008288348
## ----------End of Iteration 84 ----------
## Iteration 85 :
## Ieration 85 the differential change value of RSS (dCVR) is: 0.0007871638
## ----------End of Iteration 85 ----------
## Iteration 86 :
## Ieration 86 the differential change value of RSS (dCVR) is: 0.0007464546
## ----------End of Iteration 86 ----------
## Iteration 87 :
## Ieration 87 the differential change value of RSS (dCVR) is: 0.0007067425
## ----------End of Iteration 87 ----------
## Iteration 88 :
## Ieration 88 the differential change value of RSS (dCVR) is: 0.0006680546
## ----------End of Iteration 88 ----------
## Iteration 89 :
## Ieration 89 the differential change value of RSS (dCVR) is: 0.0006304099
## ----------End of Iteration 89 ----------
## Iteration 90 :
## Ieration 90 the differential change value of RSS (dCVR) is: 0.0005938205
## ----------End of Iteration 90 ----------
## Iteration 91 :
## Ieration 91 the differential change value of RSS (dCVR) is: 0.0005582915
## ----------End of Iteration 91 ----------
## Iteration 92 :
## Ieration 92 the differential change value of RSS (dCVR) is: 0.0005238218
## ----------End of Iteration 92 ----------
## Iteration 93 :
## Ieration 93 the differential change value of RSS (dCVR) is: 0.0004904043
## ----------End of Iteration 93 ----------
## Iteration 94 :
## Ieration 94 the differential change value of RSS (dCVR) is: 0.0004580259
## ----------End of Iteration 94 ----------
## Iteration 95 :
## Ieration 95 the differential change value of RSS (dCVR) is: 0.0004266676
## ----------End of Iteration 95 ----------
## Iteration 96 :
## Ieration 96 the differential change value of RSS (dCVR) is: 0.0003963045
## ----------End of Iteration 96 ----------
## Iteration 97 :
## Ieration 97 the differential change value of RSS (dCVR) is: 0.0003669051
## ----------End of Iteration 97 ----------
## Iteration 98 :
## Ieration 98 the differential change value of RSS (dCVR) is: 0.0003384308
## ----------End of Iteration 98 ----------
## Iteration 99 :
## Ieration 99 the differential change value of RSS (dCVR) is: 0.0003108344
## ----------End of Iteration 99 ----------
## Iteration 100 :
## Ieration 100 the differential change value of RSS (dCVR) is: 0.0002840582
## ----------End of Iteration 100 ----------
## Iteration 101 :
## Ieration 101 the differential change value of RSS (dCVR) is: 0.0002580309
## ----------End of Iteration 101 ----------
## Iteration 102 :
## Ieration 102 the differential change value of RSS (dCVR) is: 0.0002326619
## ----------End of Iteration 102 ----------
## Iteration 103 :
## Ieration 103 the differential change value of RSS (dCVR) is: 0.0002078331
## ----------End of Iteration 103 ----------
## Iteration 104 :
## Ieration 104 the differential change value of RSS (dCVR) is: 0.0001833824
## ----------End of Iteration 104 ----------
## Iteration 105 :
## Ieration 105 the differential change value of RSS (dCVR) is: 0.0001590743
## ----------End of Iteration 105 ----------
## Iteration 106 :
## Ieration 106 the differential change value of RSS (dCVR) is: 0.0001345353
## ----------End of Iteration 106 ----------
## Iteration 107 :
## Ieration 107 the differential change value of RSS (dCVR) is: 0.0001090957
## ----------End of Iteration 107 ----------
## Iteration 108 :
## Ieration 108 the differential change value of RSS (dCVR) is: 8.128309e-05
## ----------End of Iteration 108 ----------
## Iteration 109 :
## Ieration 109 the differential change value of RSS (dCVR) is: 4.60327e-05
## ----------End of Iteration 109 ----------
## Iteration 110 :
## Ieration 110 the differential change value of RSS (dCVR) is: 4.062144e-05
## ----------End of Iteration 110 ----------
## Iteration 111 :
## Ieration 111 the differential change value of RSS (dCVR) is: 6.914646e-05
## ----------End of Iteration 111 ----------
## Iteration 112 :
## Ieration 112 the differential change value of RSS (dCVR) is: 8.5722e-05
## ----------End of Iteration 112 ----------
## Iteration 113 :
## Ieration 113 the differential change value of RSS (dCVR) is: 9.70491e-05
## ----------End of Iteration 113 ----------
## Iteration 114 :
## Ieration 114 the differential change value of RSS (dCVR) is: 0.0001051343
## ----------End of Iteration 114 ----------
## Iteration 115 :
## Ieration 115 the differential change value of RSS (dCVR) is: 0.0001109407
## ----------End of Iteration 115 ----------
## Iteration 116 :
## Ieration 116 the differential change value of RSS (dCVR) is: 0.0001150325
## ----------End of Iteration 116 ----------
## Iteration 117 :
## Ieration 117 the differential change value of RSS (dCVR) is: 0.0001177807
## ----------End of Iteration 117 ----------
## Iteration 118 :
## Ieration 118 the differential change value of RSS (dCVR) is: 0.000119448
## ----------End of Iteration 118 ----------
## Iteration 119 :
## Ieration 119 the differential change value of RSS (dCVR) is: 0.0001202308
## ----------End of Iteration 119 ----------
## Iteration 120 :
## Ieration 120 the differential change value of RSS (dCVR) is: 0.000120281
## ----------End of Iteration 120 ----------
## Iteration 121 :
## Ieration 121 the differential change value of RSS (dCVR) is: 0.0001197199
## ----------End of Iteration 121 ----------
## Iteration 122 :
## Ieration 122 the differential change value of RSS (dCVR) is: 0.0001186466
## ----------End of Iteration 122 ----------
## Iteration 123 :
## Ieration 123 the differential change value of RSS (dCVR) is: 0.0001171432
## ----------End of Iteration 123 ----------
## Iteration 124 :
## Ieration 124 the differential change value of RSS (dCVR) is: 0.0001152787
## ----------End of Iteration 124 ----------
## Iteration 125 :
## Ieration 125 the differential change value of RSS (dCVR) is: 0.0001131121
## ----------End of Iteration 125 ----------
## Iteration 126 :
## Ieration 126 the differential change value of RSS (dCVR) is: 0.0001106937
## ----------End of Iteration 126 ----------
## Iteration 127 :
## Ieration 127 the differential change value of RSS (dCVR) is: 0.000108067
## ----------End of Iteration 127 ----------
## Iteration 128 :
## Ieration 128 the differential change value of RSS (dCVR) is: 0.0001052699
## ----------End of Iteration 128 ----------
## Iteration 129 :
## Ieration 129 the differential change value of RSS (dCVR) is: 0.0001023354
## ----------End of Iteration 129 ----------
## Iteration 130 :
## Ieration 130 the differential change value of RSS (dCVR) is: 9.929239e-05
## ----------End of Iteration 130 ----------
## Iteration 131 :
## Ieration 131 the differential change value of RSS (dCVR) is: 9.616616e-05
## ----------End of Iteration 131 ----------
## Iteration 132 :
## Ieration 132 the differential change value of RSS (dCVR) is: 9.297901e-05
## ----------End of Iteration 132 ----------
## Iteration 133 :
## Ieration 133 the differential change value of RSS (dCVR) is: 8.975046e-05
## ----------End of Iteration 133 ----------
## Iteration 134 :
## Ieration 134 the differential change value of RSS (dCVR) is: 8.649773e-05
## ----------End of Iteration 134 ----------
## Iteration 135 :
## Ieration 135 the differential change value of RSS (dCVR) is: 8.323595e-05
## ----------End of Iteration 135 ----------
## Iteration 136 :
## Ieration 136 the differential change value of RSS (dCVR) is: 7.997836e-05
## ----------End of Iteration 136 ----------
## Iteration 137 :
## Ieration 137 the differential change value of RSS (dCVR) is: 7.673659e-05
## ----------End of Iteration 137 ----------
## Iteration 138 :
## Ieration 138 the differential change value of RSS (dCVR) is: 7.352076e-05
## ----------End of Iteration 138 ----------
## Iteration 139 :
## Ieration 139 the differential change value of RSS (dCVR) is: 7.033971e-05
## ----------End of Iteration 139 ----------
## Iteration 140 :
## Ieration 140 the differential change value of RSS (dCVR) is: 6.720106e-05
## ----------End of Iteration 140 ----------
## Iteration 141 :
## Ieration 141 the differential change value of RSS (dCVR) is: 6.411143e-05
## ----------End of Iteration 141 ----------
## Iteration 142 :
## Ieration 142 the differential change value of RSS (dCVR) is: 6.107632e-05
## ----------End of Iteration 142 ----------
## Iteration 143 :
## Ieration 143 the differential change value of RSS (dCVR) is: 5.810054e-05
## ----------End of Iteration 143 ----------
## Iteration 144 :
## Ieration 144 the differential change value of RSS (dCVR) is: 5.518801e-05
## ----------End of Iteration 144 ----------
## Iteration 145 :
## Ieration 145 the differential change value of RSS (dCVR) is: 5.234191e-05
## ----------End of Iteration 145 ----------
## Iteration 146 :
## Ieration 146 the differential change value of RSS (dCVR) is: 4.956483e-05
## ----------End of Iteration 146 ----------
## Iteration 147 :
## Ieration 147 the differential change value of RSS (dCVR) is: 4.685875e-05
## ----------End of Iteration 147 ----------
## Iteration 148 :
## Ieration 148 the differential change value of RSS (dCVR) is: 4.422502e-05
## ----------End of Iteration 148 ----------
## Iteration 149 :
## Ieration 149 the differential change value of RSS (dCVR) is: 4.166458e-05
## ----------End of Iteration 149 ----------
## Iteration 150 :
## Ieration 150 the differential change value of RSS (dCVR) is: 3.917786e-05
## ----------End of Iteration 150 ----------
## Iteration 151 :
## Ieration 151 the differential change value of RSS (dCVR) is: 3.676487e-05
## ----------End of Iteration 151 ----------
## Iteration 152 :
## Ieration 152 the differential change value of RSS (dCVR) is: 3.442506e-05
## ----------End of Iteration 152 ----------
## Iteration 153 :
## Ieration 153 the differential change value of RSS (dCVR) is: 3.215768e-05
## ----------End of Iteration 153 ----------
## Iteration 154 :
## Ieration 154 the differential change value of RSS (dCVR) is: 2.996143e-05
## ----------End of Iteration 154 ----------
## Iteration 155 :
## Ieration 155 the differential change value of RSS (dCVR) is: 2.783453e-05
## ----------End of Iteration 155 ----------
## Iteration 156 :
## Ieration 156 the differential change value of RSS (dCVR) is: 2.577497e-05
## ----------End of Iteration 156 ----------
## Iteration 157 :
## Ieration 157 the differential change value of RSS (dCVR) is: 2.377999e-05
## ----------End of Iteration 157 ----------
## Iteration 158 :
## Ieration 158 the differential change value of RSS (dCVR) is: 2.184639e-05
## ----------End of Iteration 158 ----------
## Iteration 159 :
## Ieration 159 the differential change value of RSS (dCVR) is: 1.997021e-05
## ----------End of Iteration 159 ----------
## Iteration 160 :
## Ieration 160 the differential change value of RSS (dCVR) is: 1.814646e-05
## ----------End of Iteration 160 ----------
## Iteration 161 :
## Ieration 161 the differential change value of RSS (dCVR) is: 1.636901e-05
## ----------End of Iteration 161 ----------
## Iteration 162 :
## Ieration 162 the differential change value of RSS (dCVR) is: 1.462952e-05
## ----------End of Iteration 162 ----------
## Iteration 163 :
## Ieration 163 the differential change value of RSS (dCVR) is: 1.291706e-05
## ----------End of Iteration 163 ----------
## Iteration 164 :
## Ieration 164 the differential change value of RSS (dCVR) is: 1.121533e-05
## ----------End of Iteration 164 ----------
## Iteration 165 :
## Ieration 165 the differential change value of RSS (dCVR) is: 9.498759e-06
## ----------End of Iteration 165 ----------
mgwr.bw <- round(MGWR[[2]]$bws,1) # Nombre de voisins pour chaque prédicteur
#mgwr.bw
# Exploration des résultats statistiques
print(MGWR)## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2023-10-07 14:27:54.994629
## Call:
## gwr.multiscale(formula = prix_med ~ perc_log_vac + perc_maison +
## perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup +
## part_agri_nb_emploi + part_cadre_profintellec_nbemploi, data = data_immo_bzh_sp,
## kernel = "gaussian", adaptive = TRUE, bws0 = rep(1, 4), predictor.centered = rep(T,
## 3))
##
## Dependent (y) variable: prix_med
## Independent variables: perc_log_vac perc_maison perc_tiny_log dens_pop med_niveau_vis part_log_suroccup part_agri_nb_emploi part_cadre_profintellec_nbemploi
## Number of data points: 57
## ***********************************************************************
## * Multiscale (PSDM) GWR *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: gaussian
## Adaptive bandwidths for each coefficient(number of nearest neighbours):
## (Intercept) perc_log_vac perc_maison perc_tiny_log dens_pop
## Bandwidth 11 53 56 39 41
## med_niveau_vis part_log_suroccup part_agri_nb_emploi
## Bandwidth 56 56 56
## part_cadre_profintellec_nbemploi
## Bandwidth 56
##
## ****************Summary of GWR coefficient estimates:******************
## Min. 1st Qu. Median 3rd Qu.
## Intercept 1649.395 1723.184 1865.122 1929.827
## perc_log_vac -302.118 -298.338 -296.688 -294.292
## perc_maison -591.452 -591.063 -588.224 -583.696
## perc_tiny_log -413.370 -408.776 -402.613 -395.289
## dens_pop -219.688 -217.161 -191.024 -167.340
## med_niveau_vis 339.684 344.076 345.308 347.536
## part_log_suroccup 339.851 340.086 342.984 346.007
## part_agri_nb_emploi 98.592 102.212 103.293 104.229
## part_cadre_profintellec_nbemploi -64.223 -63.736 -59.797 -55.391
## Max.
## Intercept 1952.455
## perc_log_vac -292.470
## perc_maison -582.969
## perc_tiny_log -392.117
## dens_pop -163.384
## med_niveau_vis 353.134
## part_log_suroccup 346.272
## part_agri_nb_emploi 106.722
## part_cadre_profintellec_nbemploi -54.977
## ************************Diagnostic information*************************
## Number of data points: 57
## Effective number of parameters (2trace(S) - trace(S'S)): 14.16268
## Effective degrees of freedom (n-2trace(S) + trace(S'S)): 42.83732
## AICc value: 788.7976
## AIC value: 765.7213
## BIC value: 745.9901
## Residual sum of squares: 1837688
## R-square value: 0.8802325
## Adjusted R-square value: 0.839689
##
## ***********************************************************************
## Program stops at: 2023-10-07 14:27:58.977106
Il me semble que le nombre de plus proches voisins (number of nearest neighbours) nous indique si l’effet de notre prédicteurs relève d’un processus qui soit plus global (càd échelle du territoire) ou au contraire beaucoup plus localisé. Plus le nombre de voisins est petit plus l’effet est localisé.
Pour visualiser les Betas :
A partir de ces résultats, on peut refaire toutes les analyses et cartes réalisées avec la GWR standard.
Par exemple, les cartes des coefficients des VI par epci Breton
data_immo_bzh$agri.mgwr=MGWR$SDF$part_agri_nb_emploi
data_immo_bzh$perc_maison.mgwr <- MGWR$SDF$perc_maison
data_immo_bzh$dens_pop.mgwr=MGWR$SDF$dens_pop
data_immo_bzh$med_vie.mgwr=MGWR$SDF$med_niveau_vis
data_immo_bzh$logvac.mgwr=MGWR$SDF$perc_log_vac
data_immo_bzh$tinylog.mgwr=MGWR$SDF$perc_tiny_log
data_immo_bzh$suroccup.mgwr=MGWR$SDF$part_log_suroccup
data_immo_bzh$cadre.mgwr=MGWR$SDF$part_cadre_profintellec_nbemploi
# Réaliser la collection des cartes
par(mfrow = c(2, 4))
mf_map(x = data_immo_bzh, var = "agri.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients Agriculteurs MGWR")
mf_map(x = data_immo_bzh, var = "perc_maison.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de Maison MGWR")
mf_map(x = data_immo_bzh, var = "dens_pop.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de dens pop MGWR")
mf_map(x = data_immo_bzh, var = "med_vie.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de Médianne niveau de vie MGWR")
mf_map(x = data_immo_bzh, var = "logvac.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de Logements vacants MGWR")
mf_map(x = data_immo_bzh, var = "tinylog.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de Petits logements MGWR")
mf_map(x = data_immo_bzh, var = "suroccup.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de logement suroocupés MGWR")
mf_map(x = data_immo_bzh, var = "cadre.mgwr", type = "choro", pal= "Earth")
mf_title("Coefficients de Cadre MGWR")3.2 Régionalisation
Ici l’objectif va être de définir des sous-espaces sur la base des coefficients de la GWR, sous-espaces qui vont se caractériser par l’homogénéité des coefficients des prédicteurs au sein du sous-espace.
Thierry feuillet décrit cette méthode de cette manière :
“Ce processus de découpage de l’espace en sous-régions homogènes se nomme la régionalisation. C’est une extension de la classification classique : on y ajoute un critère de contiguité spatiale. La régionalisation est donc une classification spatiale.
Il existe plusieurs méthodes de régionalisation. Un des principes les plus répandus consiste à établir une classification à la fois sur la base de la ressemblance entre les observations, et sur leur proximité dans l’espace géographique.
Nous allons ici utiliser l’algorithme SKATER (Spatial Klustering Analysis by Tree Edge Removal), méthode proposée par Assunçao et al. (2006) et déjà appliqué dans un contexte de recherche similaire au notre par Helbich et al. (2013). Par ailleurs, une description très pédagogique de la méthode est disponible ici : http://www.jms-insee.fr/2018/S08_5_ACTE_ROUSSEZ_JMS2018.pdf
L’algorithme SKATER comporte 4 étapes (cf. doc cité ci-dessus) :
1- Constuction d’un graphe de voisinage (contiguité ou knn) 2- Pondération des liens du graphe à partir de la matrice de dissimilarité 3- Construction de l’arbre portant minimal, en retenant le lien avec le voisin le plus ressemblant pour chaque noeud 4- Elagage de l’arbre maximisant la variance inter-classes des sous-graphes
Avant de poursuivre on va donc recalculer un modèle GWR pour notre exemple breton:
## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)
best.bzh <- gwr.sel(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi,
data = data_immo_bzh,
coords = st_coordinates(st_centroid(data_immo_bzh)))## Warning: st_centroid assumes attributes are constant over geometries
## Bandwidth: 109909.9 CV score: 4736532
## Bandwidth: 177660.5 CV score: 4681812
## Bandwidth: 219532.6 CV score: 4711813
## Bandwidth: 172830.3 CV score: 4678660
## Bandwidth: 148796.9 CV score: 4668839
## Bandwidth: 136640.9 CV score: 4671999
## Bandwidth: 149754.5 CV score: 4668918
## Bandwidth: 147704.4 CV score: 4668797
## Bandwidth: 147354.2 CV score: 4668795
## Bandwidth: 147390.6 CV score: 4668795
## Bandwidth: 147393.9 CV score: 4668795
## Bandwidth: 147393.8 CV score: 4668795
## Bandwidth: 147393.8 CV score: 4668795
## Bandwidth: 147392.5 CV score: 4668795
## Bandwidth: 147393.3 CV score: 4668795
## Bandwidth: 147393.6 CV score: 4668795
## Bandwidth: 147393.7 CV score: 4668795
## Bandwidth: 147393.7 CV score: 4668795
## Bandwidth: 147393.7 CV score: 4668795
## Bandwidth: 147393.7 CV score: 4668795
## Bandwidth: 147393.7 CV score: 4668795
## Bandwidth: 147393.7 CV score: 4668795
bzh_gwr <- gwr(formula = prix_med ~ perc_log_vac + perc_maison + perc_tiny_log + dens_pop + med_niveau_vis + part_log_suroccup + part_agri_nb_emploi + part_cadre_profintellec_nbemploi,
data = data_immo_bzh,
coord = st_coordinates(st_centroid(data_immo_bzh)),
bandwidth = best.bzh,
gweight = gwr.Gauss,
hatmatrix = TRUE) ## Warning: st_centroid assumes attributes are constant over geometries
3.2.1 Première étape : préparation de la table des coefficients GWR
On fait une standardisation (centrage-réduction) pour rendre nos différents coefficients comparables.
3.2.2 Computation de l’algorithme SKATER : Définition du voisinage de chaque point
De ce que j’ai compris, Skater est un algorithme pour faire du clustering de données spatiales, en s’assurant que les clusters sont constitués d’objets contigus. Pour en savoir plus sur skater, voir ce post (skater?).
knn <- knearneigh(bzh_gwr$SDF, k = 50) # On utilise knearneigh car nous ne sommes plus sur un shape avec des polygons - on avait alors utilisé la fonction poly2nb - mais sur une matrice de points coordonnées## Warning in knearneigh(bzh_gwr$SDF, k = 50): k greater than one-third of the
## number of data points
Calibrage du coût des arêtes et de la pondération spatiale
Minimisation de l’arbre et classification
costsTree <- mstree(costsW)
plot(costsTree, coords = coordinates(bzh_gwr$SDF), col="blue", main = "Arbre portant minimal")Ici définition en 5 clusters, c’est arbitraire mais orienté par les analyses précédentes.
# Mettre spdep:: devant la fonction car le package rgeoda possède la même fonction qui fait la même chose mais pas les mêmes arguments
clus5 <- spdep::skater(edges = costsTree[,1:2], data = bzh.stand, ncuts = 5)bzhClus <- data_immo_bzh %>%
mutate(clus = as.factor(clus5$groups)) %>%
bind_cols(bzh.stand)
mf_map(x = bzhClus, var = "clus", type = "typo", pal= "Set 2")
mf_title("régionalisation Bretagne")Et grâce au code qui suit vous pouvez caractériser les clusters.
library(ggplot2)
nomVar <- c("perc_log_vac_b","perc_maison_b","perc_tiny_log_b","dens_pop_b","med_niveau_vis_b", "part_log_suroccup_b","part_agri_nb_emploi_b","part_cadre_profintellec_nbemploi_b")
clusProfile <- bzhClus[, c(nomVar, "clus")] %>%
group_by(clus) %>%
summarise_each(funs(mean)) %>%
st_drop_geometry()## Warning: `summarise_each_()` was deprecated in dplyr 0.7.0.
## ℹ Please use `across()` instead.
## ℹ The deprecated feature was likely used in the dplyr package.
## Please report the issue at <https://github.com/tidyverse/dplyr/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
##
## # Simple named list: list(mean = mean, median = median)
##
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
##
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: There were 6 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `geometry = mean(geometry)`.
## ℹ In group 1: `clus = 1`.
## Caused by warning in `mean.default()`:
## ! l'argument n'est ni numérique, ni logique : renvoi de NA
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 5 remaining warnings.
clusLong <- reshape2::melt(clusProfile, id.vars = "clus")
profilePlot <- ggplot(clusLong) +
geom_bar(aes(x = variable, y = value),
fill = "grey25",
position = "identity",
stat = "identity") +
scale_x_discrete("Effet") +
scale_y_continuous("Valeur moyenne par classe") +
facet_wrap(~ clus) +
coord_flip() +
theme(strip.background = element_rect(fill="grey25"),
strip.text = element_text(colour = "grey85", face = "bold"))
profilePlotBibliographie
Annexes
Info session
| setting | value |
|---|---|
| version | R version 4.3.1 (2023-06-16) |
| os | Ubuntu 22.04.3 LTS |
| system | x86_64, linux-gnu |
| ui | X11 |
| language | fr_FR:en |
| collate | fr_FR.UTF-8 |
| ctype | fr_FR.UTF-8 |
| tz | Europe/Paris |
| date | 2023-10-07 |
| pandoc | 3.1.1 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown) |
| package | ondiskversion | source |
|---|---|---|
| car | 3.1.2 | CRAN (R 4.3.0) |
| carData | 3.0.5 | CRAN (R 4.3.0) |
| correlation | 0.8.4 | CRAN (R 4.3.0) |
| corrplot | 0.92 | CRAN (R 4.3.0) |
| digest | 0.6.33 | CRAN (R 4.3.1) |
| dplyr | 1.1.2 | CRAN (R 4.3.0) |
| DT | 0.28 | CRAN (R 4.3.0) |
| GGally | 2.1.2 | CRAN (R 4.3.0) |
| ggplot2 | 3.4.2 | CRAN (R 4.3.0) |
| gtsummary | 1.7.1 | CRAN (R 4.3.0) |
| GWmodel | 2.2.9 | CRAN (R 4.3.0) |
| here | 1.0.1 | CRAN (R 4.3.0) |
| mapsf | 0.6.1 | CRAN (R 4.3.0) |
| maptools | 1.1.7 | CRAN (R 4.3.0) |
| Matrix | 1.6.0 | CRAN (R 4.3.1) |
| plotly | 4.10.2 | CRAN (R 4.3.0) |
| RColorBrewer | 1.1.3 | CRAN (R 4.3.0) |
| Rcpp | 1.0.11 | CRAN (R 4.3.1) |
| rgeoda | 0.0.10.2 | CRAN (R 4.3.0) |
| robustbase | 0.99.0 | CRAN (R 4.3.0) |
| sf | 1.0.14 | CRAN (R 4.3.1) |
| sp | 2.0.0 | CRAN (R 4.3.0) |
| spatialreg | 1.2.9 | CRAN (R 4.3.0) |
| spData | 2.2.2 | CRAN (R 4.3.0) |
| spdep | 1.2.8 | CRAN (R 4.3.0) |
| spgwr | 0.6.36 | CRAN (R 4.3.1) |